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Summary 


This thesis concerns exact and approximate treatments of electronic energy transfer 
in photosynthetic systems. While the methods used are completely general, their 
application is focused on the Fenna-Matthews-Olson (FMO) pigment-protein an¬ 
tenna complex, found in certain bacteria. 

The FMO complex is a trimer consisting of 24 bacteriochlorophyll (BChl) sites, each 
of which is coupled to a dissipative environment, which renders the energy transfer 
irreversible, and means that the dynamics of this transfer must be treated using 
techniques appropriate for open quantum systems. 

Recently, experimental evidence of quantum beating effects in the energy transfer 
in the FMO complex was found, and theoretical and experimental results suggested 
that quantum coherence might be observed at room temperature (300 K). One of 
the aims of this thesis is to ascertain how important coherence effects are for this 
transfer. 

Firstly, the Hierarchical Equations of Motion (HEOM) are introduced, which give 
numerically exact transfer dynamics, albeit with computational cost rising steeply 
with the size of the system. The efficient numerical implementation of the HEOM 
is discussed, and a Taylor-series integration is compared to the traditional Runge- 
Kutta method, the former proving to be more efficient. 

The accuracy of the HEOM can also be improved by representing the bath corre¬ 


lation function as a Pade series instead of the more traditional Matsubara series: 



the former converges very rapidly compared to the latter, and allows more efficient 
simulation of the dynamics. 

Results are then presented for the FMO complex, both for a 7-site subsystem of the 
monomer and for the trimer. Numerically exact results are calculated to provide 
a benchmark for cheaper approximate methods including the Redfield and Fbrster 
theories. It is found that incoherent Forster theory describes the overall features 
of the dynamics well at 300 K, suggesting that coherence has little if any effect on 
energy transfer efficiency at room temperature. 

Finally, Forster theory is used to test the effects of two phenomena on the energy 
transfer dynamics in the FMO complex: that of including vibrational structure in 
the environment (traditionally modelled as being unstructured), and that of static 
disorder due to a slowly fluctuating environment. 

It is found that energy transfer in the complex is very robust with respect to changes 
in the environment at room temperature, and that the results are largely the same 
if structure is introduced into the environment or disorder is accounted for. Rather 
than electronic coherence, it is this robustness that is advantageous for the complex’s 


biophysical role. 



Acknowledgements 


I would like to thank David Manolopoulos for supervising my Part II project 
and for his unparalleled guidance and inspiration throughout my degree. 

The members of the Manolopoulos group, Nike, Josh, Michele, Jon and Kenji, 
also helped make this year an enjoyable and productive one. Along with mem¬ 
bers of the Barford and Wilson groups, they have also provided many interesting 
conversations over many hours of coffee. 

I thank my friends for providing a distraction when I needed one, including 
Karolis, David R., Emily N., Richard, Emily B., Alan, Emma, Michael, Iain, Nat, 
Paul and Karen. 

The support of my parents has been invaluable both during my degree and 
otherwise, and I would like to thank them very much for this, along with my brother, 
Nicholas, and my girlfriend Anne, as well as my grandparents and great-aunt. 

Finally, this work is dedicated to the memory of my Nan, Joyce Wilkins, whom 


I miss very much. 



Contents 


1 Introduction [T] 

1.1 Fenna-Matthews-Olson Complex. [2] 

1.2 Open Quantum Systems. ED 

1.3 Quantum Coherence In Energy Transfer . [8] 

1.4 Summary . [9] 

2 Hierarchical Equations Of Motion [TX| 

2.1 Reduced Density Operator. [12] 

2.1.1 Superoperators And The Interaction Picture. [£3] 

2.1.2 Time-Evolution And The Partial Trace. [U] 

2.2 Influence Operator. [16] 

2.2.1 The Force Operator. [16] 

2.2.2 Evaluation Of Bath Average. [18] 

2.3 Equations Of Motion. [2U] 

2.3.1 Auxiliary Density Operators. [2T] 

2.3.2 Comments. [23] 

3 Approximate Methods [25] 

3.1 Redfield Theory. [26] 

3.1.1 Time-Dependent Redfield Theory. [26] 

3.1.2 Discussion. [2H] 

3.2 Forster Theory. [HD] 

3.2.1 Resonant Energy Transfer Rate. [30] 

3.2.2 Spectral Overlap. [33] 

4 Implementation Of The HEOM 1361 

4.1 Numerical Propagation. [36] 

4.2 Spectral Decomposition. [39] 

4.2.1 Matsubara Series. EH] 

4.2.2 Pade Series. [44] 

4.2.3 Discussion. [43] 

5 Numerical Results 1461 

5.1 7-Site Subsystem. [48] 

5.2 24-Site Trimer. [53] 

5.3 Conclusions. [HU 





























6 Applications Of Forster Theory [63] 

6.1 Structured Spectral Density. EH 

6.1.1 Results. [67] 

6.1.2 Discussion. [70] 

6.2 Static Disorder. [73] 

6.2.1 Results. [74] 

6.2.2 Discussion. EH 

6.3 Conclusions. EH] 

7 Conclusions And Further Work [T9] 

7.1 Further Work. [80] 

A Appendix [82] 

A.l The Bath Correlation Function . [82] 

A.2 The Pade Approximant. [85] 

A.3 Supplementary Information . [88] 

A.3.1 System Hamiltonians. [HE] 

A.3.2 Structured Spectral Density. [9T] 

Bibliography [92] 















Chapter 1 
Introduction 


The process of photosynthesis is one whose importance cannot be overestimated: 
photosynthetic organisms capture energy from sunlight and store it in chemical form. 
This energy is taken up by organisms higher up the food chain |TJ-3]. 

In an antenna, photons of sunlight are absorbed by pigments, generally based 
on chlorophyll molecules, and the resulting excitation is transported between these 
pigments, eventually reaching an energy “funnel”, in which the excitation moves 
over time towards pigments of lower energy, meaning that the process is irreversible. 
Its destination is the reaction centre, in which the energy is converted to chemical 
form [I]. 

This work focuses on the excitation energy transfer that occurs in the energy 
funnel, although the methods used are suitable for much more general problems. 

This transfer is extremely efficient, and it is immediately apparent that it would 
be very useful if we could design systems that were as effective at transferring energy. 
For this, it is necessary to understand the features of the system that lead to this 
efficiency. 

As an example, we will be interested in the effect of the environment surrounding 
the pigments: presumably (as with many natural systems), the environment has 
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1.1. Fenna-Matthews-Olson Complex 


been finely timed to facilitate transfer. If we were to alter this environment, what 
effect would it have on transfer rates? 


We will see in Section 1.3 that the importance of quantum-mechanical effects in 


the transfer is currently under question. By modelling the transfer using methods 
that both include and ignore these quantum effects, we can gain an idea of how vital 
they really are at physiological temperatures. 

This transfer is of current experimental and theoretical interest |3,|4l|6|-[8], and 


we now discuss a complex whose study will allow us to benchmark the techniques 
used, and to learn about the physics of energy transfer. 

1.1 Fenna-Matthews-Olson Complex 

The Fenna-Matthews-Olson (FMO) complex is a pigment-protein complex in which 
bacteriochlorophyll (BChl) pigments are bound to a protein. It is found in green 
sulfur bacteria such as C. tepidum and P. aestuarii [3j, and collects electronic exci¬ 


tation from a light-harvesting chlorosome, funnelling it towards the reaction centre. 

The complex itself is a trimer, each monomer of which was thought originally 
to comprise 7 BChl molecules [9,10, though it is now accepted that there is an 
8th, more weakly bound BChl site in addition (TTJ, [T2 . The currently accepted 
arrangement of bacteriochlorophylls is shown in Fig. 0 

The sites labelled 1, 6 and 8 are those found nearest to the chlorosome, so are 
the ones that capture the excitation energy, while those labelled 3 and 4 are closest 
to the reaction centre and lowest in energy, so it is to these sites that the energy is 
funnelled. 

There are several reasons for the popularity of the FMO complex in studies 
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1.1. Fenna-Matthews-Olson Complex 



Figure 1.1: Bacteriochlorophyll arrangement in FMO monomer. Figure created with 
Chem3D 2010, using PDB entry 3ENI [56 . 


of energy transfer: its structure is well-documented 10, 11 , and has been since 


1974 [9j, making it ideal for theoretical studies, and its solubility in water means it 
is convenient for use in experiments [4,11 . Because it has been extensively studied 


in the past, the data required to carry out our simulations is readily available 13 


Next, we turn to the problem of finding a model Hamiltonian to describe the 
excitation energy transfer, which will require an introduction to the open quantum 


system. 


3 









1.2. Open Quantum Systems 


1.2 Open Quantum Systems 

Consider a system of N “sites” (for example, BChl sites) labelled by Latin letters, 
each of which has a ground and a single excited electronic state. If |eQ refers to 
site j being in the excited state, and \gj) the ground state, we will focus on the 
single-excitation manifold, that is, states | j): 

N 

i;> = |e,>ni*>- ( L1 ) 

Using the states {|j)} as a basis, the so-called “system” Hamiltonian is given by: 

N N 

Hs = Yl huj j 1 - 7)01 + hJ M){k\ + \k)(j\). ( 1 . 2 ) 

j 

Here, hu>j is the energy of site j and fiJjk represents the coupling between two 
sites, due to dipolar interactions |2|. The eigenstates of the Hamiltonian (referred 
to as excitons) are represented by Greek letters | u), such that: 

H s \v) = huj u \v). (1.3) 


Where: 


N 


w) = ^u vj \j). 


(1.4) 


If Hs were sufficient to describe the entire process, we see that we could use a 
time-dependent state vector \ip(t)) = e lHst ^ h \ip(0)) to describe the dynamics of the 


system, and furthermore that these dynamics would be reversible (see Fig. 1.3 (a)). 


Clearly, this is not the full story, as we are not considering a completely isolated 
system. Rather, the electronic system represented by Hs is an open system, sur¬ 
rounded by an environment with which it can exchange energy. We take each site 
to be associated with its own “bath”, the vibrations of the bacteriochlorophyll and 
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1.2. Open Quantum Systems 


of the surrounding protein, and if each bath is comprised of harmonic oscillators, 
then the total bath Hamiltonian is: 

= • (i-s) 

j a ' ' 

Here, qj a is the coordinate of the a th oscillator in the j th bath, pj a is the conjugate 
momentum, m JQ the mass and c 0 j a the frequency. The system-bath interaction is: 

N N 

Hsb = Y b’) O' I 9 ^ a = %'> 0-6) 

3 a 3 

with = Ylj QjaQja the generalized force exerted on the j th site by its bath and 
Vj = \j)(j\- Before considering the effect that the environment will have on the 
dynamics of the excitation transfer, we introduce the spectral density, Jj(cj), for the 
j th site: 

J 3 M = 5^ oJ 3Q n . 6 ^ ~ W J«)- ( L7 ) 

This function has peaks at the frequencies of the bath oscillators, weighted by the 
strength of their coupling to the site j, and quantifies the system-bath interaction. 
In general, it is replaced by a smooth function of u (corresponding to an uncountably 
infinite number of oscillators). Commonly, for reasons that will be discussed later 
on in this thesis, the spectral density chosen is the Lorentz-Drude function: 




2 h Xj'yjOO 

7T 7 2 + U3 2 


( 1 . 8 ) 


Here A j = duJj(uj)/u is the reorganization energy, a measure of system-bath 

coupling strength, and 7 j is the characteristic frequency of the bath, which satisfies 
dJj (ca) 


du 


= 7 j) = 0. Fig. 


1.2 


illustrates this spectral density. 


The harmonic oscillators of the bath cause a kind of quantum-mechanical friction, 


one effect of which is exactly the same as classical friction: energy is dissipated from 
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1.2. Open Quantum Systems 



Figure 1.2: Lorentz-Drude spectral density of Eqn. 
35 cm -1 , 7 j = 106.1 cm -1 . 


( 1 . 8 ). 


Parameters used: A j = 


the electronic system, giving irreversible energy transfer. 

There is, in addition, the purely quantum-mechanical effect of decoherence. The 
vibrational states associated with each electronic state may be excited or de-excited 
(we refer to them as absorbing or emitting a phonon), which alters the site energy 
|37|. This in turn causes the coherence between two sites to decay, as the phase e ld 
of one state | k) with respect to another, | j), is altered when the energy of a site 
changes, and when these phases are randomized, their average value is zero. Over 
time, quantum-mechanical oscillations will be damped. 


Both of these phenomena are illustrated in Fig. 1.3, which shows the dynamics 


of a closed system and of a system with exactly the same Hamiltonian Hs, but 
coupled to a bath of harmonic oscillators. 

Since a larger reorganization energy A j gives a stronger interaction between bath 


and system 27 , we will expect stronger dissipation and more decoherence. For the 


same reason, higher temperatures also lead to greater decoherence. 
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1.2. Open Quantum Systems 


(a) (b) 



Figure 1.3: Comparison of the dynamics of (a) a closed and (b) an open quantum system, 
showing the effects of dissipation and decoherence. 


One further property of the open quantum system must be mentioned: that of 
Markovianity. Due to the statistical nature of the environment, the force acting on 
a given site is a stochastic function of time. 

If this stochastic process were Markovian, then at any given time no memory 
of its past behaviour would be required to predict its behaviour in the future; only 
its current value would be needed. There is no reason to assume a priori that 
the behaviour of the function is Markovian, so that an accurate description of the 
dynamics should allow for non-Markovian effects. 

The techniques used to deal with open quantum systems are applicable not only 
to biological systems such as the FMO complex, but also to many others, such as 


electron and proton transfer 25-27 , as well as damping of Rabi oscillations due to 


excitation in quantum dots 29 . The theoretical tools described in this work are 


thus very versatile, and useful in many fields of current scientific interest. 
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1.3. Quantum Coherence In Energy Transfer 


1.3 Quantum Coherence In Energy Transfer 

We might intuitively expect that for the FMO complex at a physiological tempera¬ 
ture of around 300 K, there would be a very strong interaction between environment 
and system (due to thermal excitation of bath oscillators), leading to very fast de¬ 
coherence. 

Fairly recently Engel et al. j4] observed quantum beating effects in a 2-dimensional 
electronic spectroscopy study of energy transfer in the FMO complex at a cryogenic 
temperature of 77 K, and Collini et al. j5| observed similar effects at an ambient 
temperature of 294 K for a system similar to the FMO complex, leading to the 
suggestion that quantum coherence is important for energy transfer, for example in 
tunnelling through energy barriers |6|. 

This result was followed by a number of theoretical studies by Ishizaki and Flem¬ 
ing [6 20,38 with the aim of accurately modelling the energy transfer dynamics in 
the FMO complex. 

Importantly, using the Hierarchical Equations of Motion (which will be discussed 
in this thesis), the dynamics were simulated at both 77 K and 300 K for a 7-site 
subsystem of the monomer |6|. At both temperatures oscillations were observed, 
suggesting that even at room temperature, quantum coherence was seen in the 
energy transfer, thus sparking renewed interest in the FMO complex. 

One of our aims in the following work will be to appraise whether or not quantum 
coherent effects are necessary for efficient energy transfer at room temperature. 


In the literature, it is conventional to discuss excitation transfer dynamics by 
showing time-dependent site populations |6| (i.e., in the site basis (|j)}), and we 







1.4. Summary 


will use this convention in the work to follow. 

A number of timescales are important in the study of the FMO complex: elec¬ 
tronic coherence lasts for a number of femtoseconds (around 500 fs at 77 K, and 
around 200 fs at 300 K), so that a numerical simulation up to 1 ps will capture the 
effects of coherence, if appropriate. 

After the coherence has decayed, the population dynamics show a simple decay, 
like those predicted by a rate equation, and eventually the populations reach a 
steady state. This occurs at around 10 ps at 77 K, and around 5 ps at 300 K, and 
will be observed in a numerical simulation up to 15 ps. 

A number of experiments have been carried out which show that the fluorescence 


lifetime of bacteriochlorophylls in the FMO complex is on the order of 1 ns 14 . This 


means that we need not consider loss of excitation energy via spontaneous emission 
on the timescales used in our simulations, as such events are suitably rare. 


1.4 Summary 

A brief survey of the remainder of this thesis is now appropriate. Firstly, in Chapter 
[2] the Hierarchical Equations Of Motion (HEOM) for calculation of exact energy 
transfer dynamics will be introduced, and the advantages and disadvantages of using 
these equations will be described. 

Next, in Chapter [3] we will introduce several cheaper methods that can be used 
to obtain approximate results for the dynamics, the Redheld and Forster theories. 

Chapter [4] contains a description of some methods used to implement the HEOM 
efficiently, while doing so more accurately than in previous work. 
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1.4. Summary 


Chapter [5] is dedicated to numerical results for the 7-site and the 24-site FMO 
complex (the latter of which has not had any exact results published yet) from the 
HEOM and from the approximate methods described in Chapter |4j It is found that 
at 300 K, the incoherent Forster theory describes the general features of the energy 
transfer quite well. 

Chapter [6] then covers applications of Forster theory: we find the effects of using 
a structured spectral density instead of the Lorentz-Drude function, and of static 
disorder due to slow fluctuations in the protein environment. 

Finally, Chapter [7] concludes. 
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Chapter 2 

Hierarchical Equations Of Motion 


There are several methods for finding a numerically exact description of the 
dynamics of an open quantum system, many of which have been in use for some 


decades. The Feynman-Vernon influence functional 17 is based on the path integral 


formulation of quantum mechanics 18 , wherein the dynamics of the system alone 


are considered, having integrated out those of the bath. 


A more recent method is the self-consistent hybrid method 27 , in which the 


spectral density is not considered as a continuum, but is discretized, giving a number 
of “bath modes”. The dynamics of a number of these modes is treated exactly along 
with the system, while the dynamics of the rest of the modes is treated classically: 
the number of modes treated exactly is increased until the results converge. 

In this Chapter we present the Hierarchical Equations of Motion (HEOM), which 


were formulated originally by Tanimura and Kubo 15,16 using the influence func 


tional as a starting point. We work instead with operators 20 , deriving the influ¬ 


ence operator from a consideration of the statistical properties of a bath of harmonic 
oscillators. 


Section 2.1 introduces the reduced density operator, and gives the general equa¬ 


tion of motion for this operator. In Section 2.2, the influence operator is derived, 
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2.1. Reduced Density Operator 


and in Section 2.3, a hierarchy of auxiliary density operators is introduced, giving 


exact equations of motion. 

2.1 Reduced Density Operator 

An ensemble of quantum systems is most appropriately described not with wave- 
functions but rather, using the density operator formalism. If H is the Hamiltonian 
for a system and {|hq)} the eigenstates of H, then the density operator is given by: 


p = 


( 2 , 1 ) 


where pj is the probability that a system in this ensemble will be found in the state 
\4>j). Particularly important is the situation in which the system is a member of 
a canonical ensemble, so that pj = e" 13 ^ /Zp, where f3 = 1/ksT, €j is the energy 
eigenvalue corresponding to eigenstate \4>j), and Zp = tr[e~ BH ] is the canonical 
partition function. In this case, the canonical density operator is pp = IZp. 

A diagonal matrix element of the density operator, ((pk\p\4>k), gives the popula¬ 


tion of state |0*;) 19 , and by cyclic permutation of operators within a trace, the 
average value of an operator A is given by: 


(A) = = tr[Ap\. 

j,k 


( 2 . 2 ) 


For a system described by a single state vector, if this state is not an eigenstate of 
the Hamiltonian then it will evolve in time, and similarly, from the Time-Dependent 
Schrodinger equation, if a density matrix does not commute with the Hamiltonian, 


we will observe dynamics 19 , given by ihdp(t)/dt = [H,p{t)). 
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2.1. Reduced Density Operator 


2.1.1 Superoperators And The Interaction Picture 


It proves useful now to introduce some definitions. Firstly, the Liouvillian superop¬ 


erator C is defined such that: 


CO = [. H,0}. 


(2.3) 


Further, we can define two superoperators A ~ 4 and A^~ such that A^O = AO 
and A^O = OA. Then the commutation and anticommutation superoperators are 


given by [21]: 


i x = A* - A<~ => A x O = [A, o\, 


A° = + A*~ => A°o = {d, 6}. 


(2.4) 


(2.5) 


A useful relation is the following, which can be verified by differentiating both 
its left and right hand sides: 


e iCt/hQ — e iHt/hQ e ~iHt/h 


( 2 . 6 ) 


In order to derive an equation of motion for the density operator, we first split 
up the total Hamiltonian as H = H 0 + V (and thus, the Liouvillian as C = Co + Cy)- 
Then, in the interaction picture with respect to Hq, denoted by a tilde following the 


notation of 20 , the state vectors and operators are given by (pp. 172-173 of 1221): 


1 4>(t)) = e iflot/h \ <j,{t)) = e iAot/n e~ i6t/n \(j)(0 )), 


(2.7a) 


0(t) = e idot/n O. 


(2.7b) 


Inserting Eqn. (2.7a) into the time-dependent Schrodinger equation gives: 




( 2 . 8 ) 
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2.1. Reduced Density Operator 


and this gives the expression for the time-derivative of the density operator in the 
interaction picture: 

j t m = ( 2 . 9 ) 

Solution of this differential equation will give an expression for the time-dependence 
of the density operator for the system described by Hamiltonian H. 


2.1.2 Time-Evolution And The Partial Trace 


Solution of (2.9) is not as straightforward as it first appears. To understand this, 


we first integrate the equation to give a recurrence relation: 


p(t) = m ~j l J dti£v(h)p(h). 


( 2 . 10 ) 


The density matrix at a given time is determined by an integral over density 


matrices in the past, and can be expressed as an infinite series by expanding (2.10): 


pit) = m+ 


k =1 


k rt ftk -i 

-n) L*'-L 


dtk£v(ti) ■ ■ ■ £v(tk) 


P(0). (2.11) 


There is a notable similarity between this and the so-called perturbation expan¬ 


sion 


18 . By analogy, we can interpret (2.11) as a sum over alternative possibilities. 


A general term in the sum represents the effect of Hsb on the system at times 
ti,t 2 , ■ ■ ■ ,tfc- Liouvillians at two different times do not necessarily commute, and so 
the order in which they are applied is important. 

We introduce the time-ordering operator T■ For a product of two operators, the 
action of T is: 


T0i(ti)0 2 (t 2 ) — Oi(ti)O 2 (t 2 )0(ti — t 2 ) + O 2 (t 2 )Oi(h)0(t 2 — ti). (2.12) 
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2.1. Reduced Density Operator 


Here, O (t) is the Heaviside step function (equal to 1 when t > 0 and 0 when 
t < 0). When T acts on a product of n operators, there will be a total of n\ terms on 
the right hand side (this is the number of ways that the operators could be ordered). 


Applying T to (2.11) so that the influences of the Cy{t) are experienced in 


chronological order (thus preserving causality) gives: 


Pit) = p / (o) + rJ^ 


k =1 


— ^ J dt\ ... J dtkCv{ti) ... Cy{tk ) 

pt 


M o) 


I «»)■ 


This is the power series for the exponential function, so can be rewritten: 


(2.13) 


p(t) = T exp ( -- / dhCv(ti) ) p(0). 


(2.14) 


The above discussion is entirely general, and is suitable for any quantum-mechanical 
system, including the system-plus-environment supersystem we are considering. 
However, the environment has an infinite number of degrees of freedom, making 


it impossible to implement (2.14) in the form shown 


A popular method of dealing with this problem is to take a partial trace; the 
reduced density operator is a trace, over all bath degrees of freedom, of the total 
density operator: 

Ps(t) = trs[p(*)]- (2.15) 


Finally, we assume that the initial density matrix is given by p( 0) = Ps{0)Pb,/3, 
where Pb,p — e~^ HB /Zb is the canonical density matrix for the bath (note that no 
tilde is used in this case, since ps,p = Pb,p and the canonical density operator is 
unchanged in the interaction representation). 
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2.2. Influence Operator 


This is justified in an electronic excitation process by invoking the Franck- 


Condon principle 20 : when electronic excitation occurs, the baths remain in their 


equilibrium state because nuclear motion is slow on the timescale of electronic reor¬ 
ganization. 

With the definition = tr B[0{t)pB,p\, the time-evolution of the reduced 

system density matrix in the interaction picture is given by: 


Ps(t) = T ( exp -- / ) Ps(0). 


(2.16) 


This expression is extremely important, and will be revisited in the work of 
Chapter [3j 


2.2 Influence Operator 


Evaluating the bath average in (2.16) will give an expression that depends not on the 


bath operators but only on terms that give the influence of the bath on the system. 
In order to achieve this, some observations must be made about the bath. Firstly, 
we set Hq = Hs + Hb and V = Hsb , so that in the remainder of this Chapter the 
interaction Liouvillian is Cy{t) = CsB(t). 


2.2.1 The Force Operator 

The operator £j(f) = gj a qj a (t ) gives the force exerted on the j th site by the bath 
oscillators a associated with this site. It is our aim here to find an explicit form for 
£j(t), and to show that the force exerted on a site by its bath is a random, Gaussian 
force, which will allow the bath average to be performed simply. 
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2.2. Influence Operator 


In the interaction picture, differentiating Eqn. (2.7b) gives 


fo {t ) = iA, o(t). 


(2.17) 


Using Hq = H s +Hb gives dq } Jt)/dt. = and dqqjtj/dt = 


whose solution is: 


pj a ( 0) 

Qja[t) = Qja{ 0) COS [0Jj a t) H--- Sill (LUj a t) . 


17lj a UJj a 


(2.18) 


From this we infer 23 that the force operator at time t is determined by the 


positions and momenta of the bath oscillators at time 0, and so we now show that 
these positions and momenta are Gaussian variables. That is: 


(f(Qj a ))p oc 



dq ja f{qja)e K< b 2 “, 


(2.19) 


for any function ) of qj a in the Schrodinger representation, and that a similar 
relation will be true for functions of momentum. 

Recalling the definition of the bath average, the trace can be evaluated in the 
position basis. The trace is then given by in, a dqj a (qja | • • • | qjo)- Since the canon¬ 
ical density matrix is a product pb ,/3 = Tij Pj,P = integration over each 

coordinate qy a t gives the partition function Zy a t, except for coordinate qj a : 


( f{qja))p — / dqj a (qj a \f(qj a )pj a y\qj a )/Zj C 


— ry / dqj a f(qj a )(qj a \pj at /3\qj a ). 


( 2 . 20 ) 


J ja J—oo 


The canonical density matrix element ( qja\Pja,p\qja) is perhaps most easily eval- 


(t), 
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2.2. Influence Operator 


uated by the path integral formulation of quantum mechanics 18 , which give^J 


(f(Qj a ))p = yr [ ^Qjoif isijot) exp f- m ^ a tanh () q 2 


J ja J —oo 


h 


ja 


( 2 . 21 ) 


The positions and momenta of thermal bath oscillators are stochastic variables 
(due to the statistical nature of the problem) with Gaussian distributions. Since a 

, we 


linear combination of variables with Gaussian distributions is also Gaussian 24 


have the result that the force operator £j(t) acting on site j has such a distribution. 


2.2.2 Evaluation Of Bath Average 


We can now eliminate the dependency of Eqn. (2.16) on the operators of the bath. 
Using the identity ( AB) X = \{A X , B °} + j{A°, B x } + \ [A, B} x for the commutator 


of a product of operators 21 gives: 


N 




N 


E(&w°^w x +&w x ^w o )- 


( 2 . 22 ) 


We can infer that the interaction Liouvillian also has a Gaussian distribution 


with respect to a bath average, as it will also be a linear combination of the force 

operators. This distribution means that the statistical properties of CsB(t) are fully 
1 This result is obtained by noting that the canonical density matrix e~P H is analogous to the 
propagator for “thermal time” t = —ifiTi, and treating it using the path integral method 

for evaluating such propagators. 
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2.2. Influence Operator 


determined by its autocorrelation function 24 




vAtr v j ( S y+ 

wr(h(tri(s)°) fl vA‘‘)*+ 

vAtr^AtriiW^/iW 0 ) 


(2.23) 


Noting that, for example, ^£ 7 (t) 0 £j(s)°y = tr (t)°£j{s)°p-j.p] , this trace will 
have four terms, and can be simplified (by cyclic permutation of the operators 
within the trace) to give 2 £j(s)}^ • The remaining three bath averages can 

be carried out (with the third and fourth vanishing) to give: 


^SB(t)CsB(s)^ — ... 


N 

£V)(t) x (a rj -(t-s)y;-(s) x + ia ij -(t*s)^(s)°) . (2.24) 


Here a r j(t) = \ &(0)}y and ia itj (t) = ] , so that 


oij{t) = Oi r j(t) + ) is equal to ^£j(t)£j(0)y , the autocorrelation function of 

the force exerted on the system by the bath. 

Using the identity for a Gaussian process x(t) [24], 


exp ( / du/(u)z(u; 


= exp 


dfi / dt 2 f(ti)f{t 2 ) (x(ti)a;(t 2 )) , 


(2.25) 


19 






2.3. Equations Of Motion 


we finally obtain 120 [ 


N / 1 rt pt i 

Ps(t) — T exp ( ^ / dt\ I dt 2 V j (t i y a rJ (ti - t 2 )Vj{t 2 )' 


'o ./o 


+ zaij(ti - t 2 )Vj(t 2 ) 0 )ps(0). (2.26) 


The operator acting on ps(0) is the influence operator |20|. In the position 
representation, this becomes the influence functional of Feynman and Vernon flT,28|, 


which expresses the time-propagation of the reduced density matrix in terms only 
of statistical quantities of environment variables. 


It is shown in Appendix A.l that the bath correlation function a-j(t) can be 
expressed in terms of the spectral density, such that Jj(-uj) = —Jj(co): 


Oij[t) — / Jj(w) [coth(/37ku/2) cos(u;t) — i sin(u;t)] duj. 

Jo 


(2.27) 


Thus, the influence operator encodes the effects of the baths (through the spectral 
densities) and of the temperature (through the coth.(/3hu/2) term) on the system. 


2.3 Equations Of Motion 

With an explicit form for the bath correlation function we can now find a 

differential equation to evolve the reduced density matrix. We take: 


K 


a jit) = p j0 5+(t) + ^Pjk e ljkt 


(2.28) 


k =1 


Here, S + (t) is a one-sided delta function, which has all of the properties of the 
delta function (pp. 58-60 of 1221), but the integration limits in its definition are 


Note that Eqn. (|2.25 ) has been adapted here due to the time-ordering; the operator T allows 


us to change the integration limits and removes the factor of 
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2.3. Equations Of Motion 


(0,oo) instead of (— 00 , 00 ). For example (and most importantly), J 0 °° 8 + (t)dt = 1. 
We also use the notation pj k = a,j k + ibj k , where a,j k and bjk are real. 

The reasoning behind this expression is given in Chapter [4| and the advantage 
of this sum of exponentials will quickly become apparent in the following work. 


2.3.1 Auxiliary Density Operators 


Using (2.28), we can write the time-evolution of the density operator as 16 


Ps(t) = T JJ exp (Jp J q dU0j(U)%o(U)^ x ... 

JJexp J dti J dt 2 $j(t 1 )e~^ tl ~ t2) 6 jik (t 2 )^ As(0), (2.29) 


where we have introduced the notation, 




Oj, k (t) = ia jk Vj(t) x - b jk Vj{t)°. 


(2.30a) 


(2.30b) 


We define a set of auxiliary density operators (ADOs) a n (t), indexed by the 
matrix n whose elements are rij k : 


<7n(t) = T]^[exp ^^2 J dU0j(ti)0yo(U)^ X ... 

( 1 pt \ n jk 

x... 

exp (Jp f f *20 J Vi)e- 7 « ( “- bl # At (i2)Us(O). (2.31) 
The ADO with n = 0 is the reduced density operator (RDO). The reason for 


these ADOs being introduced is made clear when the time-derivative of Eqn. (2.31) 
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2.3. Equations Of Motion 


is taken 20 : 


n(0 — 'y ] o{t) n jk r YjkJ &n{t) + 


qrCTnlt) = 

at 


j,k 


Y + n jk 9 jtk {t)a Ujk _{t )J . (2.32) 

j,k 

The ADOs <7 n (i) and cf lljfc± (t) have matrices n that are identical, except that for 
the latter the element rijk is replaced by rijk ± 1. For each ADO, we dehne: 


N K 

At n = njk- 

j k 


(2.33) 


And now, noting that dO(t)/dt = e lCot ^ n dO/dt + jC Q e lCot ^ h O , the time-derivative 


(2.32) can be rearranged and rewritten in the Schrodinger picture: 


^d n (t) = -J^sd n (t) - Y “ ib joVj*Vj + WjfcTjfc) ^n(t) + • • • 

j,k 

Y (^ X(3 n Jfc+ (t) + in jk p jk Vjd n . k _ (t) + in jk p* k d njk _ (t)V 3 ^ . (2.34) 

j,k 

We see now that these ADOs form a hierarchy: a given ADO’s level in the 
hierarchy is given by its value of -M n , with only one operator (the reduced density 
operator) in the lowest level, N x K operators in the first level, and so on. Each 
ADO is coupled to operators on the levels above and below. 

An interesting qualitative picture of the HEOM is given by reference to creation 
and annihilation operators (pp. 136-139 of 1221): to each element of n we can assign 
a “bath mode” labelled by jk. Then, a given ADO <r n has rijk quanta in the bath 
mode jk. 

It is then possible to see that a given level of the hierarchy corresponds to a 
certain total number of quanta, and that each ADO is connected to the ADOs that 
can be formed from it by creating or annihilating one quantum. 
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2.3. Equations Of Motion 


By propagating these ADOs through time (we see in 2.3.2 that a finite number 
of these operators can be used), the time-evolution of the diagonal elements of the 
RDO (the site populations) can be found. 


2.3.2 Comments 

The remainder of this Chapter is given to a discussion of several important points 
about the HEOM. 


• We have made no approximation in this derivation, so have a method that 
is numerically exact and able to describe the effects of quantum-mechanical 
coherence in energy transfer and of non-Markovian system-bath interactions. 

• For all ADOs but the RDO, the initial condition is that all elements of these 
matrices be equal to zero. 

• Other than the RDO, the ADOs cannot be system density matrices, as their 
traces are not conserved. 


• In principle, the hierarchy continues to infinity. However, this is not possible 


in practice. Integrating (2.34): 


cr. 


L (i) = [ dt x • e -(^s+E i ,4^o^y/-^oC/v/+n jfc7 , fc ))(t-t 1 ) x 


E (A^.+M + in jk p jk Vja njk _ (ti) + in jk p* k a njk _ . (2.35) 
j,k 

For ADOs above some level A f of the hierarchy, the values of rij k will be so 
large that exp (— ^ krtjkljkit — H)) will decay very rapidly and thus will be 


proportional to 5 + (t — 1 1 ). This gives the result 16,20 that for these ADOs, 
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2.3. Equations Of Motion 


the time derivative does not depend on ADOs in any other levels, so their 
elements will remain at their initial value of zero. 


• We obtain converged results by carrying out calculations with an increasing 
number of levels AT, until increasing this number has no further effect on the 
resulting dynamics. 


• For a calculation using J\f levels of the hierarchy, the number of ADOs required 



(JV+ (NAT))! 
A/1(N A')! 


(2.36) 


For some calculations, this number can be extremely large, meaning that the 
HEOM can be very computationally expensive, requiring a lot of time and 
memory to run. We will wish to implement the equations as efficiently as 
possible; this will be discussed further in Chapter |4j 
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Chapter 3 

Approximate Methods 


Although the HEOM provide a very powerful method for calculating energy 
transfer dynamics, the large computational cost incurred means that we will need 
to seek an approximate method with a much lower cost, but with dynamics that are 
as accurate as possible, since we will be unable to use the HEOM for several of the 
physical features of the FMO complex that we wish to investigate. 

There are a great number of methods that are suitable candidates, both historical 
and modern. Some of these, such as the Forster (1948) |32| and Redheld (1957) |30 
theories, treat some physical quantity as a perturbation; these two methods will be 
presented in this Chapter. 

Other methods include the Zeroth Order Functional Expansion quantum master 
equation (ZOFE, 2011) |34l, in which a key operator is treated as being independent 


of the bath, and which has recently been applied both to the 7-site FMO system 46 


and the 24-site trimer 47 


This more modern method will not be considered in this thesis, as preliminary 
calculations (not reported here) showed the ZOFE to be less accurate than simple 
Forster theory for the FMO complex at 300 K. 
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3.1. Redfield Theory 


The existence of so many methods, all with different conditions for applicability, 
means that it is necessary for us to decide which we will use in our investigations, 
and which we need not consider any further. We turn now to the task of introducing 
two techniques whose validity we will test in Chapter [5} 


3.1 Redfield Theory 


The perturbative nature of both the Redfield and Fdrster theories 35 means that 


the validity of each depends on some parameter being small enough to justify it being 
treated as a perturbation. In the case of the Redfield theory, the coupling between 
the system and bath is assumed to be weak, characterized by reorganization energies 
A j that are small (compared to dipolar couplings J 3 k)- Thus, HgB is taken as the 
perturbation to the Hamiltonian. 

In qualitative terms, we might consider an electronic system evolving under 
its own Hamiltonian Hg, which is then weakly coupled to an environment. The 


phonons destroy phase information 37 , assisting the relaxation of the system to an 


equilibrium state. 


3.1.1 Time-Dependent Redfield Theory 


We take as a starting point Eqn. (2.16), with Cy(t) = CsB(t)- Taking the time- 
derivative gives: 


^Ps(t) = (^-jT£ S B(t)exp ^-1 I dtiC S B(ti)Jj ps( 0 ) 


—| A J 

~ — * (£sB(t)^^ps(0) — j2 J dti (3.1) 


h 


with the second line being the series expansion up to second order in the perturbation 
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3.1. Redfield Theory 


Hsb■ The hrst term on the right vanishes, since it is an average of a Gaussian 


variable, which is zero 24 . 


The correlation function appearing in the integrand is given explicitly by (2.24). 


We also replace ps(0) on the right hand side by ps(t), assuming that the weak 
perturbation causes only a negligible change in ps(t) between time 0 and time t. 


Then, expanding the commutators and anticommutators 20 
d 


dt 


Ps(t ) = Y / dt^ajit - ti)V j (t)V j {t±)p s (t) + 


- aj(t - tjVjitJpsitfVjit) — a*(t — t 1 )V j (t)p s (t)V j (t 1 )+ 


a*(t - ti)p s (t)Vj(t 




• (3-2) 


We choose to evaluate the density matrix in the exciton basis, where Hs\p) = 
Tico^p), because this means that the matrix elements of Vj(t) have the form: 


(p\e i6st/n Ve- i6st/n \u) = e^\p\V\v)e 


—iuj u t 


(3.3) 


The density matrix can be transformed between exciton and site representation 


by using the matrix U of Eqn. (1.4) that diagonalizes Hs 


37 


In finding the matrix elements of (3.2), the completeness relation is used: , | p!) (p'\ 

1. For the hrst term on the right-hand side, the matrix elements are given by: 

1 X 

- 72 X] / dtiUjit - tiWVMVjitdPsWW) = 
n j Jo 


l 

¥ 


yvW dt\ • ctjit - x 

Jo 

(p | Vj | k) (k| Vj\p) (p | p s (t) | v). (3.4) 


3 K M 


Dehning a “partial Fourier transform”, 


aj[u;t\ = / aj(ti)e liJtl dti, 

Jo 


(3.5) 
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3.1. Redfield Theory 


and changing the variable of integration, Eqn. (3.4) can be rewritten as: 


e i(.Un-Wv)t 

—^ 2 — 2^ 2^ 2^ s ™' 2^t\(t j '\j)(3\ K )( K \j)W)pi*',A t )- (3.6) 


Here, pyy(t) = {p!\p(t)\v'), and (j\n) = U Kj , as in Eqn. (1.4) 


Treating the other three terms in (3.2) in the same manner gives the Time- 


Dependent Redfield equation 37,38 : 


^Pn,u{t) ^ ] Ry,,v,n'y(t)Pn'y(f)■ (3-7) 


fl' V' 


The key quantity in the Redfield theory is the matrix R^vyyit), which governs 
the relaxation of the electronic system to equilibrium. It is given by: 


Rfj,,v,y y (t) — T ul ^y(t) + r^/ ^ v y (t) 


8fv' ^H,K,K,n' {t) 8 flfl 1 ^V,K,K,v’{ 0 - ( 3 - 8 ) 


Here T^yyit) is the time-dependent damping matrix, and its elements are 38,39 : 


r w,yy(t) = ■^^{p 1 \j){j\v){p!\j){j\u')a j [u vl -ay;*]. (3.9) 


The time-dependence of the relaxation matrix (3.8) highlights the non-Markovian 


nature of the Time-Dependent Redfield theory: Eqn. (3.5) is responsible for memory 


effects 36,39 , and using Eqn. (2.28), this function is given by 39 : 


aqlo;; t\ = 


X 1 — e (*w-7 jfc )i 
;*] =Pjo + 2^Pjk- 


k=l 


7j k 


(3.10) 


3.1.2 Discussion 


The following points are salient when considering Redfield theory: 





















3.1. Redfield Theory 


Due to the assumption of weak system-bath coupling in Eqn. (3.1), the full 


dynamics are not captured: the perturbation Hsb is proportional to a linear 
combination of bath coordinates, so can only induce single-phonon transitions 
in the bath [35] (there will only be non-zero matrix elements between bath 
states whose energy levels are adjacent). 


In its original form 30 , the Redfield theory was Markovian. This theory 


can be derived simply from the Time-Dependent Redfield theory by assuming 


that e“ 7 * ~ 7 _1 h + (t) in Eqn. (3.5): qualitatively, by time t, the integrand has 
already decayed to zero, so that the upper bound of the integral can be set to 
infinity without appreciating its value. 

• The only change made to the Time-Dependent theory to derive the Markovian 


theory is to replace ctj[u);t] by oij[u>] = aq[u;; t —y oo] 36 , giving a relaxation 
matrix that is time-independent. 

When implementing the Time-Dependent Redfield theory to evolve the density 


matrix through time, the relaxation matrix (3.8) must be calculated at each 


timestep, whereas it is not dependent on time in the Markovian Redfield case 
so only needs to be calculated once. 


• More recently, a modified Redfield theory has been developed 35,40 , which 


can interpolate between the Redfield and Forster theories. However, this the¬ 
ory ignores off-diagonal density matrix elements in the exciton basis, so cannot 
be used to describe dynamics in the site basis [3j. 
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3.2. Forster Theory 


3.2 Forster Theory 


Forster theory, or Forster-Dexter theory 32,33,41 , was originally formulated to 


describe resonant energy transfer (RET) between two electronic states (a donor and 
an acceptor) |42|. 


More recently, Jang et al. introduced a generalized Forster theory 43 , which 


takes into account situations in which environmental phonons have not relaxed to 
equilibrium before energy transfer occurs. 

In this Section, we derive the original Forster theory, and then explain how to 


calculate the resulting energy transfer rate constants. 


3.2.1 Resonant Energy Transfer Rate 

For clarity, we begin by considering two sites, labelled 1 and 2, each of whose 
electronic states are associated with a continuum of vibrational modes, so that when 
there is no interaction between the sites, the vibronic states are, 

W(Ei,E 2 )), (3.11) 

where Ej is the vibronic energy of site j. Hq is the unperturbed Hamiltonian: 

2 

Ho — hujj\j)(j\ + Hb + H S b, (3-12) 

3 =1 

such that E 2 )) is an eigenstate of H 0 : 

{^(Ei, E2)\Ho\' l P(Ei, E 2 )) —E 1 + E 2 . (3.13) 

Now, we add an interaction J (which has no diagonal matrix elements between 
vibronic states) to the Hamiltonian (note that this gives Ho + J = H, the full 
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3.2. Forster Theory 


Hamiltonian): 

j= £«*( |j><*| + |fc>0D- (3.14) 

By construction, this interaction is weak (its matrix elements are smaller than the 
reorganization energies A j) and treating it as a perturbation, it will induce vibronic 
transitions, leading to energy transfer between sites. 

We take site 1 to be the “donor” and site 2 the “acceptor”, with the former being 
electronically excited initially and the latter being unexcited, giving an initial state 
denoted e 2 )) (where e\ is a vibrational level of the excited electronic state of 

site 1 and likewise e 2 for the ground electronic state of site 2 ). 

Invoking Fermi’s Golden Rule ( |3l], p. 178 of 122]) gives the probability of being 
in some state \i/j(Ei, E 2 )) as a function of time. Since this probability is the square 
modulus of a single probability amplitude, we are ignoring phase information and 
any theory built on this framework will be incoherent. The probability is given by: 

P(E U E 2 , t; 61 , € 2 , 0) = |Me,, e 2 )\j\^(E h g 2 ))| 2 ■ - , (3 . 15 ) 

where A E = + e 2 ) — (Ei + E 2 ) is the energy difference between the initial and 
final vibronic states. The total rate of transfer from site 1 to site 2 is given by 
integrating over all final energies E\ and E 2 (where E\ is a vibrational level of the 
ground electronic state of site 1 and likewise E 2 for the excited electronic state of 
site 2 ): 

P\—>2 (t‘i GT2; 0) = 

dE t jT dE 2 m (l , £2 )| J\i,(E u £ 2 )>| 2 ■ 4Sm2 (fJ t/2ft) . (3.16) 

The integrand contains a function of the form sin 2 (hixt)/x 2 (where k is a con¬ 
stant), which in the limit as t —y 00 is proportional to the delta-function. We thus 
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3.2. Forster Theory 


make the approximation 31 : 


4sin 2 (AF't/2ft) 

AE ~ 2 


7 Tt ( A E 


2nt r , . x 

= iA (AB) - 


(3.17) 


and use the property (p. 60 of [22]): 


5(AE) = S(e i + e 2 — E 1 — E 2 ) = / d(Uu) ■ 5{Uu + e 2 - E 2 )5(hu — ei + E x ) 

J —oo 

/ oo 

du ■ 5(E 2 — [e 2 + hoj])S(Ei — [ei — hu>]). 

-OO 

(3.18) 

This gives: 


/ oo /»oo /»oo 

dca / dE t / d£ 2 |(?/i(ei,e 2 )| j)?/i(£i,£ 2 ))| 2 x 

-oo «/ 0 «/ 0 


h(£ 2 - [e 2 + - [ei - M), (3-19) 


and we see, since A E = 0 (due to the approximation in Eqn. (3.2.1[)), that energy 


is conserved, and from Eqn. (3.19), that Hut is the energy transferred from donor to 


acceptor | 31| . 

Finally, the rate constant for energy transfer is found by using the definition 
F f( x )d( x — a)dx = f(a ) and taking the time-derivative of the probability: 


fcf^ 2 ( e iA 2 ) = 27T / du\('ip(e 1 ,e 2 )\j\ip(e 1 -huj,e 2 + huj))\ 2 . (3.20) 


Since the probabilities depend linearly on time, the rate constants are time- 
independent. Next, we will rewrite the integrand in terms of the more familiar 
functions of open quantum systems. 
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3.2. Forster Theory 


3.2.2 Spectral Overlap 


We now assume that the system is initially at equilibrium, so that we do not precisely 
know the initial vibronic energy; rather, we know that the probability density of 
being in state |-0(ei,e 2 )) is given by g(ei)g(e 2 ), where g(e) = e~ ,3e /Z and Z = 
J 0 °° e~P e de is the canonical partition function. The total rate of transfer is given by: 


2 = 2?r 


doj 


de i / de 2 ■ g(e 1 )g(e 2 )x 


|(V’( e i, e 2 )\H\il>(ei - huj,e 2 + huj))\ 2 . (3.21) 


Within the Born-Oppenheimer approximation, the vibronic state is given as a 
product of electronic and vibrational states: 


mE u E 2 )) = \^)\xi(E 1 ))\x2{E 2 )), (3.22) 


where | Xj(Ej)) is the vibrational state of site j. We also take |1) as the initial 
electronic state, with only the donor excited, and |2) as the final electronic state, 
with only the acceptor excited. Then, using (1| J|2) = hJ\ 2 . 


k^ 2 = 2nh 2 \J 12 \ 2 I du [ j dei • 5f(ei)|(xi(ei)|xi(ei - ^,o;))| 2 J x 

/ /*oo \ 

de 2 ■ g{e 2 )\{x 2 (^ 2 )\X 2 i ^2 + kiu))\ 2 \ . (3.23) 

The fluorescence and absorption spectral line shapes, Fi[u;] and A 2 [uj\ respec¬ 
tively, are defined in terms of Franck-Condon factors S 2 (ej, ej±hco) = \(xj{ej)\Xj( e j^ z 

M) I 2 

coo 

Fi[uj] = 2nTi / de\ ■ g(ei)S 2 (ei,e\ — Tilu), (3.24a) 

Jo 

COO 

A 2 \lS\ = 2^1% / de 2 ■ g( y e 2 )S 2 (e 2 , e 2 + hu>). (3.24b) 

Jo 
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3.2. Forster Theory 


The rate constant is given by: 


i r°° 

k (_+ 2 = 7H-M 2 / F^A^ujjduj. 

J —oo 


(3.25) 


All that is needed now is a method to calculate these spectra using known pa¬ 
rameters such as the bath correlation functions oij(t ), site energies huij and reorga¬ 
nization energies hXj. 


This connection is made in |35,44, 55 , where the fluorescence and absorption 


spectra are given by Fi[oj] = Fi(t)e lut dt and A 2 [uj\ = A 2 (t)e lU}t dt. In the 
time-domain, the fluorescence and absorption line-shape functions are evaluated 


using the cumulant expansion method 44,45 , giving: 


Fi(t) = e "®( ul - Al ) i "9l( t ) j 

A 2 {t) = g-^+Ai )t-92(t)' 


(3.26a) 


(3.26b) 


The function gj(t) is defined: 


9j(t) = ~2 [ dt\ [ dt 2 ■ oij(t 2 ). 
h Jo Jo 


(3.26c) 


h\j is the reorganization energy of site j, the energy difference between the 
excited vibrational level reached by a Franck-Condon transition and the ground 


vibrational level, in a given electronic state 36 


Parseval’s theorem states that the overlap integral in Eqn. (3.25) can be carried 


out in either domain, so that with some rearrangement, and assuming that energy 
transfer in a donor/acceptor pair is not affected by the presence of other sites, the 
Forster rate constant for transfer between two electronic states \j) and | k) is given 
by the equation (3ft denotes the real part) (35): 


kf^ k = 21 «/,-fc| 2 3ft 


)jk\ 


Fj(t)A k (t)dt 


(3.27) 
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3.2. Forster Theory 


Fig. 3.1 illustrates the spectra Fi[o;] and A 2 [u;], showing their overlap at both 
77 K and 300 K. The effect of the f un ction gj(t ) in F\(t) and A 2 (t) is to broaden 
the spectra, from the delta-functions observed if ctj(t) = 0 (that is, in the absence 
of baths). 


(a) 77 K 


(b) 300 K 




Figure 3.1: Overlap of the fluorescence spectrum of site 1 and the absorption spectrum 
of site 2 in the FMO complex at (a) 77 K and (b) 300 K. 


At the higher temperature, the bath correlation function is larger, the spectral 
lines are broader and there is more overlap between them, thus increasing the rate 
of transfer from site 1 to site 2. 

Chapter [5] will test the applicability of the Forster theory by comparing the pre¬ 
dicted dynamics to the exact dynamics of the HEOM. The quality of the predictions 
of Forster theory will allow us to determine whether or not coherence is important 
in the physics of the FMO. 
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Chapter 4 

Implementation Of The HEOM 


The number of ADOs required for a HEOM calculation, Eqn. (2.36), can be very 


large, and the computational cost of propagating these matrices through time is high. 
The first Section of this Chapter explains how the HEOM are solved numerically, 
with an eye towards maximizing efficiency. 

The second Section explains why the bath correlation function a(t) can be rep¬ 
resented as a sum of exponentials (which is essential for the HEOM), and explores 
both the Matsubara and the Pade spectral decompositions, highlighting the greater 
accuracy that can be achieved by using the latter. 


4.1 Numerical Propagation 


Eqn. (2.34) provides a set of coupled differential equations for all of the ADOs. If 


a(t) is a vector containing all of the elements of all ADOs, the equations can be 


rewritten: 


j t o{t) = A-2.it). 


(4.1) 


The second form follows from the fact that the differential equations do not 
contain the elements of the ADOs to any order greater than the first: A is a linear 
operator. We now look at the numerical solution of these equations. 
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4.1. Numerical Propagation 


Originally, the well-established fourth-order Runge-Kutta (RK4) method was 
chosen for numerical integration (Sec. 15.1 of |48j). If timestep St is used for the 
integration, then the error is of the order of St 5 . However, an alternative numerical 


solution was noted. The formal solution to Eqn. (4.1) is a(t + St) = exp (ASt)a(t) 
which can be Taylor-series expanded: 


M 


St 


a(t + St) = V — r A m ff(t) + 0(St M+1 ). 
m\ — 


(4.2) 


m =0 


If M — 4 is chosen, the error will again be on the order of St 5 . The Taylor series 
method has an advantage that is not immediately obvious: it leads to a reduction 
in the memory required to run a simulation. 

The reason for this is that the RK4 method involves calculation of four vectors 
k x , k 2 , k 3 and k 4 . It is possible to use only two of these vectors, with one used for 
permanent storage of k x and the other updated so that it will variously store the 
remaining three vectors. 

On the other hand, the Taylor-series method requires that only one such vector, 
k x , is used. The program will loop through four steps, each time calculating A A ■ k x 
and adding this to the current a(t). 

The RK4 method requires that the program for implementing the HEOM stores 
one more vector that contains as many elements as there are in all of the AD Os 
(that is, J\f x N 2 , where JV is as in Eqn. (2.36)). This can lead to a significant 
reduction in memory usage. 

The time-derivative of each ADO depends on the values of operators on the levels 
above and below its own. The number of ADOs in a given level is potentially huge, 
so that an efficient indexing system is vital: before the time-evolution begins, an 
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4.1. Numerical Propagation 


indexing matrix is set np, each of whose rows corresponds to a single ADO, and 
whose columns contain references to the operators to which the ADO of this row is 
coupled, for ease of use at runtime. 

For a given problem, we will have a system Hamiltonian H$ and a bath character¬ 
ized by a(t), a sum of exponentials. By increasing the number of these exponentials, 
onr a(t) will approach the true function. This number and the number of hierarchi¬ 
cal levels may be increased until the dynamics converge and are numerically exact. 

Within the program to carry out the HEOM, both the indexing routine and 
the calculation of the time-derivatives were parallelized: running on multiple CPUs 
allowed a faster, more efficient program. 

The routine calculating the time-derivatives involves (when a(t) contains no 
(5 + -function) the calculation of 6 matrix products for each ADO, or 6N 3 scalar 
multiplications. However, for the commutator and anticommutator of V 3 with a 
general operator O, a general matrix element is given by: 

{m\V^O\n) = S mj (n\0\n) - (m\0\j)S jn , (4.3a) 

(m\V°0\n) = S mj (j\0\n) + (m\0\j)8 jn . (4.3b) 

Using this formula instead of explicitly calculating the commutators and anti- 
commutators will give only 2N 3 scalar multiplications (in finding the commutator 
£b n (f)), and led to a substantial speeding up of the program, thus increasing effi¬ 
ciency. 
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4.2. Spectral Decomposition 


4.2 Spectral Decomposition 


The Lorentz-Drude spectral density is frequently used for calculations because it 
allows the bath correlation function to be written analytically as a sum of decaying 
exponentials, meaning that the HEOM can be used to benchmark our calculations. 


Traditionally, the integral in Eqn. (2.27) (or equivalently (A.10)) is carried out 


as described below to give the Matsubara series [6,16,25.49 ; this series converges 
very slowly to an exact result, and so we describe an alternative, the Pade series, 
which is found to converge very much faster. 


4.2.1 Matsubara Series 


Inserting the definition of J{oS) from Eqn. ( 1 . 2 ), with subscripts dropped for clarity, 


into (2.27) and making use of the fact that the integrand is an even function gives 


the integral: 


a(t) = 


q f\h 


dco 


7r 


—iut 


co coth(/3Tico/2)e lut toe 

to 2 + y 2 + co 2 + y 2 


(4.4) 


This can be solved straightforwardly using contour integration: the integrand has 
poles at c o = ±* 7 , and (for the first term) at sinh(/3/kn/2) = 0 =»■ co = ±i2irj//3h, 
with j an integer. As we are only interested in a(t) for t > 0 , we use the contour 
below, with R —» 00 : 
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4.2. Spectral Decomposition 



Re (a;) 


The integral over the semi-circle vanishes (pp. 113-115 of |50 ), so that the 
integral over the real line is given by a sum of residues: 


(t) = —2yA hi < Res + Res 

• u=—i r ) ‘ ^ u=—i2'Kk/ fih 


k= 1 


oj coth( flhuj/2)e 


-iut ue -^t 


oj 2 + y 2 


+ 


oj 2 + y 2 


= 7 A h | cot 




— 1 > e 




E 


47A u k 


0 -»kt 


rr P K - V) 


(4.5) 


Here, v k = 2nk/f3Ti is called a Matsubara frequency and we have a sum of 
decaying exponentials, as desired, it is possible to rewrite cot( 0 ) as an infinite sum 
by applying the residue theorem once again (pp. 131-133 of [501), giving the form 
used by Ishizaki and Fleming | 6 |: 


E 


2 7 2 


a(t) = ^:i 

P 1 


:P h "l \ ,-Tt __ - x y , 


/? 


7/2 


-1 


fe=l 


^ -7" 


(4.6) 


It is from this expression that a truncation is suggested: a certain number K 
of exponential terms are retained, and for k > K, a Markovian approximation is 
used: it is assumed that u k is large enough that u k e~ Vkt ~ S + (t), or equivalently that 
these terms decay quickly enough that they are essentially zero for t > 0 , and only 
contribute at t = 0 . 
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4.2. Spectral Decomposition 


In practice, one would successively increase K until converged results were found 


for the dynamics. Fig. 4.1 shows aft) for some values of K. 



Figure 4.1: Convergence of the Matsubara series at 77 K with varying number of expo¬ 
nential terms K, for the Lorentz-Drude spectral density. Inset: Magnified bath correlation 
functions at short times. 


We see that a large value of K would be required to fully converge the HEOM 
calculations, and so we next turn our attention to a better method of representing 
aft). 


4.2.2 Pade Series 


Faster convergence can be obtained by finding the Pade approximant to the Bose- 


Einstein function, fsose{ x ) = (1 — e ~ x )^ 1 51 . 52 1 . That is, we find a rational function 
that approximates fBosefx)- We will focus on the [N-l/N] approximant (for reasons 
that will be discussed at the end of this Section), which uses a fraction of the form: 


Eg ^ 

Zlohx* 


(4.7) 
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4.2. Spectral Decomposition 


In Appendix |A.2, it is shown that the Bose-Einstein function can be approxi¬ 


mated 52 


N 


fBoseipc) ~ T + 2x ^ ^ 

nr / ' 


Vj 


where the r]j can be found using: 


V j =[N^ + -N 


3 , r \ iL v ,'(a-cj 


=i s . 


2 ’ 
3 


(4.8) 


4 = 1,2,...,1V 


(4.9) 


2 J nh,K - ej) ’ 

Here, £.,■ = 2/cj and £j = 2/c.,-, where are the positive eigenvalues of matrix 
A and Cj the positive eigenvalues of matrix A. A is a 2N x 2N matrix, and A a 


2N — 1 x 2N — 1 matrix, with elements 52 


A rn:n 


-^m,nd=l 


V / (2m + l)(2n + l) 


A-m/n. 


-'m,nd=l 


\J (2m + 3)(2n + 3) 


(4.10) 


There will be N parameters and IV — 1 parameters Cj- 


Inserting (4.8) into (A.10) gives, with Uj = C,j//3h (hereafter, the Pade frequen¬ 


cies): 


, x 2 7 Ah , ue~ iujt 
a(t) = —— / du- 


1 1 2u rjj 


7T 


1 Z.UA V—■V 

+ 7+«*£ 


( w 2 + 7 2 ) [Phu ' 2 ' phj^uj 2 + v] 

This integral can be evaluated using the same contour used to derive the Mat- 


(4.11) 


subara series, to give, by analogy with Eqn. (4.6): 


2A 




2 t 7 j- 7 2 _ n//3h \ t 2A ^ 27^77 t 


N1 ^-7 2 ' 2 


£ A V ^ 


(4.12) 


Fig. 4.2 shows how the Pade series for a(t ) changes as the N of [. N-l/N] is 
varied. The improved convergence is very noticeable: while the Matsubara series 
required a large number of terms to give a good approximation, the Pade series 


requires a much smaller number. 
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4.2. Spectral Decomposition 



Figure 4.2: Convergence of the [N-l/N] Pade series at 77 K with varying N, for the 


Lorentz-Drude spectral density. Inset: Magnified bath correlation functions at short times. 


The code used to produce this figure was supplied by David Manolopoulos. 


Using the Pade decomposition, we can hope to predict the dynamics of energy 
transfer with an accuracy that might not be afforded to us by the Matsubara de¬ 
composition: HEOM simulations can be run with increasing number of exponentials 
until convergence is achieved. 

On the other hand, because the Matsubara series requires so many more expo¬ 
nential terms to achieve convergence, it is likely that at some point, the computa¬ 
tional expense will prohibit more exponentials from being used, before convergence 
is achieved. 


4.2.3 Discussion 

The Matsubara and Pade series are not the only methods of representing a{t ) as 
a sum of exponentials: there are several alternatives, including the Meier-Tannor 
decomposition [53j, in which the spectral density is fitted to an expression of the 
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4.2. Spectral Decomposition 


form: 


■/(-) = ?£ 


u 


7“i ((cu + Dfc ) 2 + r|)((o; - Dfc ) 2 + r 2 ) 


(4.13) 


Another possibility is to numerically integrate Eqn. (2.27), and then to fit the 


numerical data to an expression of the form (2.28). It may be appropriate, in 


this case, to allow both the prefactors and the frequencies in this expression to be 
complex. 

Both of these methods involve numerical fitting, but since the Pade series al¬ 
lows efficient analytical convergence to the exact correlation function (and efficient 
numerical convergence to the exact energy transfer dynamics) for a Lorentz-Drude 
spectral density, it is this method that we use. 

As well as the [N-l/N] approximant, the [N/N] and [A/’+l/A/’] approximants 
were also considered in the original literature |52|. The latter two require that 


a (5 + -function approximation be introduced, whereas the former involves no such 
approximation. 

We also briefly consider the high-temperature limit in order to make the connec¬ 
tion with literature in this held (6 16 . The imaginary part of the bath correlation 
function is independent of temperature, cq(t) = — Ayhe -7 *. However, for the real 
part, using lim^ 0 coth(/37iu;/2) = 


2Ay f 00 cos (cat) 2A , 
lim a r (t) = —f / doj ; = — e~^. 
n/3 J 0 cu 2 + y 2 0 


(4.14) 


This high-temperature expression is purely classical, and we might consider fur¬ 
ther terms (due to the Pade or Matsubara expansion) as “quantum corrections” to 


this 16 


Interestingly, even at high temperatures cq(f) contains a factor of h, so is inher- 
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4.2. Spectral Decomposition 


ently quantum-mechanical. This term can be expressed as the Fourier transform of 
J(oo), so can be related to quantum dissipation, whereas the /^-dependence of a r (t) 
relates it to equilibrium fluctuations, which are classical at high temperatures. 
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Chapter 5 

Numerical Results 


We have now laid the theoretical foundations necessary to carry out calculations 
for the 7-site FMO system and the 24-site trimer, at a temperature of 300 K. Calcu¬ 
lations at a temperature of 77 K have also been carried out, but are not shown here 
due to the length constraint for this thesis, as well as the fact that this will allow 
us to focus on results at physiological temperature, which are more relevant to this 
work. 


We will use the same parameters as in the literature |6|, so that h\ n = 35 cm -1 
and h'y n = 106.1 cm -1 for each site. These values are found by fluorescent Stokes 
shift experiments [6,55 . 


For the 7-site FMO, the system Hamiltonian used is that of Adolphs and Renger 


13 , given in Appendix A.3, Eqn. (A.21). The couplings hJjk between sites were 


found using the transition dipole moments of these sites and assuming the protein 
to provide a dielectric medium, while the site energies hujj were calculated using the 


interaction between the sites and charged amino acid residues 13 


For the 24-site FMO, we follow the theoretical study of Ritschel et al. 47 , whose 


ZOFE quantum master equation calculations used two different sets of site energies. 


One of these was found by Schmidt am Busch et al. 12 using a method similar to 
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4.2. Spectral Decomposition 


that of Adolphs and Renger 13 , and the other by Olbrich et al. 56 using molecular 


dynamics simulations and electronic structure calculations. The Hamiltonian is 


given in Appendix A.3, Eqn. (A.22), with Eqn. (A.23) giving the intra-monomer 


couplings, (A.24) giving the site energies for the two different cases 12,56 and 


(A.25) the inter-monomer couplings. 


The initial conditions we use will reflect the physics of the complex: bacteri- 
ochlorophylls 1, 6 and 8 are closest to the chlorosome, and it is these that will be 
most likely to be excited initially 12j 47 . Thus, for the 7-site system we will carry 
out simulations with excitation beginning on sites 1 or on 6, and for the trimer, on 
sites 1, 6 or 8. 


The results are colour-coded according to the key in Fig. |5.1[ with the population 
of each site denoted by the colours shown. 


— BChl 1 
BChl 2 

— BChl 3 
BChl 4 
BChl 5 
BChl 6 
BChl 7 

— BChl 8 


Figure 5.1: Key for numerical results. 

We set these results out as follows: firstly, those for the 7-site system are pre¬ 
sented, both for 1 ps and 15 ps (i.e., on the timescale of the electronic coherence 
and on the timescale of equilibration), and then those for the trimer. 

These results will then allow conclusions to be drawn about which approximate 
method is most appropriate for our further work. 
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5.1. 7-Site Subsystem 


5.1 7-Site Subsystem 


Figs, T2 and A3 compare the dynamics of Time-Dependent Redfield and Forster 
theories to those of the HEOM at 300 K, up to 1 ps, with the initial excitation on 


site 1 and on site 6 respectively. Figs. 5.4 and 5.5 show the same dynamics up to 


15 ps. In each case, the populations of only four sites are shown, as the rest of the 
populations stay close to zero. 

At 300 K, the HEOM dynamics fully converged with 2 Pade exponential terms 
in a(t) and 4 levels of the hierarchy. At lower temperatures, the weaker system- 
bath interaction means that fewer levels of the hierarchy are required, but more 
exponential terms 

In these results, electronic coherence is observed on short timescales. The for¬ 
mation of a coherent superposition between two states | j) and | k) is expected when 
the magnitude of .7^, the coupling between them, is comparable to |o ij — 04 . |, their 


energy gap. Consulting Eqn. (A.21), it is for this reason that coherence is observed 


between states |1) and 12), as well as between | 6 ) and |5). 

Whether the first or the sixth site is initially excited, the steady state is the 
same, but if site 6 is initially excited then this state is reached more quickly. Of the 
7 sites, BChl 6 has the highest energy, and we might thus intuitively expect that 
transfer away from this site is faster, which is borne out by the results. 

Calculations were also carried out with the Markovian Redfield theory, and com¬ 
pared to the Time-Dependent Redfield results, the only difference was that the 
short-time oscillations were less pronounced using the Markovian theory. After the 
envelope of these oscillations decayed, both theories were quantitatively identical. 
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5.1. 7-Site Subsystem 


(a) Time-Dependent Redfield 



Time / ps 


Figure 5.2: Energy transfer dynamics for the 7-site FMO subsystem at 300 K, with initial 
excitation on site 1. In each case, the solid lines show the HEOM result and the dashed 
lines show (a) Time-Dependent Redfield, (b) Forster results. 
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5.1. 7-Site Subsystem 


(a) Time-Dependent Redfield 




Figure 5.3: Energy transfer dynamics for the 7-site FMO subsystem at 300 K, with initial 
excitation on site 6. In each case, the solid lines show the HEOM result and the dashed 
lines show (a) Time-Dependent Redfield, (b) Forster results. 
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5.1. 7-Site Subsystem 


(a) Time-Dependent Redfield 




Figure 5.4: Long-time energy transfer dynamics for the 7-site FMO system at 300 K, 
with initial excitation on site 1. In each case, the solid lines show the HEOM result and 
the dashed lines show (a) Time-Dependent Redfield, (b) Forster results. 
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5.1. 7-Site Subsystem 


(a) Time-Dependent Redfield 




Figure 5.5: Long-time energy transfer dynamics for the 7-site FMO system at 300 K, 
with initial excitation on site 6. In each case, the solid lines show the HEOM result and 
the dashed lines show (a) Time-Dependent Redfield, (b) Forster results. 
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5.2. 24-Site Trimer 


Forster theory’s performance tends to be worse in general than that of Time- 
Dependent Redheld theory, especially at shorter timescales, as is particularly no¬ 


ticeable in Fig. 5.3 


Both Redheld and Forster theories predict the equilibrium populations very well 


at 300 K. However, as seen in Fig. R5 the perturbative methods predict a “bump” 
in the dynamics that is not seen in the HEOM simulation. 

Overall, the reasonable accuracy of Forster theory’s predictions at 300 K suggests 
that even if coherent effects are important at low temperatures, they do not seem to 
be very much so at room temperature: energy barriers can be surmounted thermally 
and therefore tunnelling contributions have little effect. 


5.2 24-Site Trimer 

In this Section, we present accurate results for the full FMO trimer up to 1 ps, as well 
as results up to 15 ps that are converged with respect to number of hierarchy levels, 
but not fully converged with respect to number of exponentials in a(t) (however, it 
is known that these results are very accurate). These results have not previously 
been published. 

The reason that the 15 ps results are not converged with respect to number of 
exponentials is one of computation: up to the shorter time, converged results using 
the HEOM required a large amount of time and memory, and a prohibitively larger 
time would be required to observe the fully converged steady states. 

However, by using the same bath correlation function for all methods, we can 
carry out a rigorous comparison, which will be illuminating. 
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5.2. 24-Site Trimer 


Figs. 5.6 and 5.7 compare the Time-Dependent Redfield and the HEOM results, 
using respectively the site energies of Olbrich et al. (Olb) and of Schmidt am Busch 


et al. (SaB), while Figs. 5.8 and 5.9 repeat this comparison for the Forster and the 
HEOM results. 


Fig. |5.10 compares the results of the HEOM up to 15 ps to those of the Time- 


Dependent Redfield, and Fig. 5.11 compares the results of the HEOM to those 
of Forster theory on the same timescale. Both of these Figures use the Olb site 
energies. 

When BChl 8 is initially excited, the decay is largely exponential, with no co¬ 
herence effects, due to the fact that population transferred to BChl 1 is rapidly 
transferred onwards to BChl 2, and sufficient population for a noticeable coherent 


superposition is not built up 47 


Considering the Time-Dependent Redfield simulations, we can see that this 
method performs quite well for both sets of site energies, capturing the features 
of the transfer dynamics including the oscillations. 

At 300 K, Forster theory is quite accurate overall, and in cases where it is not as 
much so, the differences are generally over-predictions of transfer rates rather than 
under-predictions. The results are somewhat worse than those of Redfield theory: 
aside from its neglect of the oscillations (whose importance we consider in the next 
Section), the agreement with exact dynamics tends to be poorer for Forster theory. 
Both theories nevertheless predict the long-time dynamics well, as was seen for the 
7-site model. 
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5.2. 24-Site Trimer 


(a) Initial Excitation On Site 1 





Figure 5.6: Energy transfer dynamics for the FMO trimer at 300 K using Olb site energies, 
with HEOM (solid lines) compared to Time-Dependent Redfield (dashed lines). Initial 
excitation on (a) site 1, (b) site 6, (c) site 8. 
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5.2. 24-Site Trimer 


(a) Initial Excitation On Site 1 





Figure 5.7: Energy transfer dynamics for the FMO trimer at 300 K using SaB site 
energies, with HEOM (solid lines) compared to Time-Dependent Redfield (dashed lines). 
Initial excitation on (a) site 1, (b) site 6, (c) site 8. 
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5.2. 24-Site Trimer 


(a) Initial Excitation On Site 1 





Figure 5.8: Energy transfer dynamics for the FMO trimer at 300 K using Olb site energies, 
with HEOM (solid lines) compared to Forster theory (dashed lines). Initial excitation on 
(a) site 1, (b) site 6, (c) site 8. 


57 
































5.2. 24-Site Trimer 


(a) Initial Excitation On Site 1 





Figure 5.9: Energy transfer dynamics for the FMO trimer at 300 K using SaB site energies, 
with HEOM (solid lines) compared to Forster theory (dashed lines). Initial excitation on 
(a) site 1, (b) site 6, (c) site 8. 
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5.2. 24-Site Trimer 


(a) Initial Excitation On Site 1 




Figure 5.10: Comparison of long-time dynamics predicted by the HEOM (solid lines) and 
Time-Dependent Redfield (dashed lines) theories for (a) initial excitation on site 1, (b) 
initial excitation on site 2. The Olb site energies are used. 
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5.2. 24-Site Trimer 


(a) Initial Excitation On Site 1 




Figure 5.11: Comparison of long-time dynamics predicted by the HEOM (solid lines) and 
Forster (dashed lines) theories for (a) initial excitation on site 1, (b) initial excitation on 
site 2. The Olb site energies are used. 
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5.3. Conclusions 


5.3 Conclusions 

In this Chapter we have presented numerically exact electronic energy transfer 
dynamics for both the 7-site and the 24-site FMO models, and have used them 
to benchmark approximate calculations using the Time-Dependent Redfield and 
Forster theories, in order to identify the method most suitable for further use. 

Comparing the Forster and Redfield theories, although the Forster theory pre¬ 
dicted the general shape of the dynamics and tended to give reasonably good quan¬ 
titative agreement, the predictions of the Redfield theory were more accurate. Both 
methods, however, predicted the long-time decays observed, and the steady state, 
well. 

The fluorescence lifetime of the bacteriochlorophyll sites is on the order of nanosec¬ 
onds, whereas, as we have seen in this Chapter, steady states are reached in a 
timescale on the order of tens of picoseconds. Whether or not electronic coherence 
effects are included in the model of the energy transfer, this steady state is reached 
in the same amount of time. 

The fact that Forster theory predicted the dynamics at 300 K so well, in terms of 
the important features (transfer rates, equilibrium populations and the timescale at 
which a steady state is reached) is extremely interesting, and allows us to conclude 
that although quantum coherence effects may be present in the energy transfer at 
physiological temperature, they appear to be unnecessary for efficient transfer. 

With this conclusion, we are prepared to further investigate the physics of the 
FMO complex, using the Forster theory, which has been shown to give fairly good 
results at temperatures of biophysical interest and is computationally extremely 
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5.3. Conclusions 


cheap compared to all other methods described in this thesis, including the Time- 
Dependent Redheld theory. This will be the subject of the following Chapter. 
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Chapter 6 

Applications Of Forster Theory 


There are a number of physical features of the FMO complex which cannot be 
investigated using the HEOM, due to the computational demand. Forster theory, 
on the other hand, requires very little in terms of computational resources, so can 
be used for these investigations with ease. 

The popularity of the Lorentz-Drude spectral density utilized thus far, within 
the modelling of open quantum systems, is due in large part to the fact that it allows 
the bath correlation function to be written as a sum of exponentials. 

In reality, there is no reason to expect that this should give us physically cor¬ 
rect results: in a pigment-protein complex, the featureless spectral density we have 
been using does not describe the vibrational environment of a bacteriochlorophyll 
pigment: for example, vibrations due to certain bonds (for example, O-H) might be 
expected to couple strongly to a site. 

A biological system is generally exquisitely tailored, with its physical parame¬ 
ters striking a balance which, if disturbed past their tolerance, could have drastic 
effects on its functionality. The spectral density contains all information about the 
environment, and so altering this function is equivalent to altering the environment 
that interacts with the bacteriochlorophyll. 
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6.1. Structured Spectral Density 


In Section 6.1, we explore the effect of using structured spectral densities, and in 


doing so we hope to learn how robust the energy transfer dynamics are to a change 
in environment. 

The site energies given for the FMO complex are in fact averages over some 
distribution: due to the slow fluctuation of the protein environment, the actual set 
of site energies for a given complex will differ from this average, a phenomenon 


known as static disorder 57 . 


The effect of this disorder, simulated by drawing the site energies from a Gaussian 
distribution, on the energy transfer dynamics of the complex will be the subject of 
Section 16.21 

In this Chapter, we carry out simulations only at a physiological temperature of 
300 K, as this is more biologically meaningful, and it is at this temperature that we 
have shown the Forster theory to give reasonably accurate results for the important 
features of the energy transfer dynamics. 

6.1 Structured Spectral Density 

The molecular dynamics simulations used by Olbrich et al. to find the site energies 


of the FMO trimer were also used to calculate a spectral density for this complex 56 


58 . Fig. 6.1 shows an average spectral density for the 7-site subsystem, normalized 


to give a reorganization energy of A = 35 cm , in agreement with experiment. 

The simulations carried out gave the real part of the bath correlation function 


in the form: 


K 


a r ,j(t) = Cjk cos(ui jt t)e : ' !l , 


( 6 . 1 ) 


k= 1 


where K = 15 and the values of the Cjk, tOjk and 'fjk are given in Appendix A.3 
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6.1. Structured Spectral Density- 



Figure 6.1: Structured spectral density for 7-site FMO complex, using parameters given 


m 


58 


Eqn. ( |A.26 ). From this, Fourier inversion of the real part of (2.27) gives the spectral 
density^ Jj(oj) = | tanh(/3huj/2) f 0 °° a r j(t) cos(ut)dt, or: 

tanh(/3huj/2) K 




£ 


Cjk^fjk 


+ 


C'jkHjk 


* ^ \'ijk +( w - + ("+%) 


(6.2) 


In order to use this spectral density in our calculations, we need an expression 
for the full bath correlation function, so must find 


Oii,j(t) — — / dojJj(u)sva(ut). 

Jo 


(6.3) 


Using tanh(j3huj/2) = j— — 1+e 1 /3 ^ gives an integrand that is an odd function, 
so that the integral can be rewritten: 


OL. 


hi 


l _i — 


C-jk'Yjk 


+ 


Cjk'Jjk 


7T 


°° duj . sin M) 

k=1 ^-oo 1 + 2 jk +{ui - uijk) 2 Ijk + (u + tijkY 

UU r°° ( uj 2 + u) 2 fc + 7 2 fc ) sin (ut)f Fer mi{fJhuj) 

2_s Cjkljk / dui r _ 2 - / - ^ ir ., 2 , c. , \2) • ( 6 - 4 ) 

k=l 


hjk + ( w ~ hjk + ( w + %*0 2 ] 


factor of ^ is included to give agreement with j58j, although the scaling to give the correct 
A means that this factor is unimportant 
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6.1. Structured Spectral Density 


Here, fFermi( x ) — (1 + e x ) 1 is the Fermi-Dirac function. As with the Bose- 


Einstein function in 4.2.2, there is an [N-l/N] Pade approximant that can be used 
for f Fe 


[x 


f Fermi \ x ) 


1 


N 


-2xJ2 


Vi 


i=i 


x 2 +e' 


For the Fermi-Dirac case, the matrices A and A are given by: 


(6.5) 


_ ^m,nd=l 

mn / . , . = 

\J (2m — l)(2n — 1) 


^ _ ^m,nd=l 

\J (2m + l)(2n + 1) 


( 6 . 6 ) 


Using the Pade expression gives the following integral, where A denotes the 
imaginary part: 


^i. j (t) 


, K N 

—yy 

fBU z —^ z —^ 

^ k=l 1=1 


CjkHjkTU^S 


doj-e~ luJt x 


v 3 + tf k u + 


bfjk + ( w _ %fc) 2 ] [t jk + ( w + Ujk) 2 ) [ix 2 + V 2 


. (6.7) 


This integral can be carried out using the usual semicircular contour. The result 
is, with Q jk = iu jk - - y jk : 


K N 


<*«(«) = 4 EE' CltmSilt 


v? - 

' k=i i=i \ 1 


ft? 


'jk 


CjkVl^jk Q* t 

—I— --- P jk 

- (%) 2 

2^ j kC j Pli{p j k\- - vf) Ul t 

h-n%\ 2 


. ( 6 . 8 ) 


The real part can also be written in terms of exponentials: 


K 


a 


k=l 


(6.9) 


Since each term in ( 6 . 2 ) gives rise to a peak in the spectral density, by excluding 


certain terms we can now find out whether removing the corresponding peaks has 


an effect on the energy transfer dynamics. 
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6.1. Structured Spectral Density 


6.1.1 Results 


In order to inform the cases we investigate in this Section, we will wish to know 


the electronic transition frequencies between sites. Fig. 6.2 shows these frequencies 


along with the structured spectral density. 



Figure 6.2: Inter-site transition energies for the 7-site FMO complex. 


The investigation of Olbrich et al. using the ZOFE |46| led to the conclusion 
that high-energy vibrational modes do not have any effect on the energy transfer 
dynamics, since the energy of these modes is far from being in resonance with any 
transitions. 


The spectral density in 6.1 shows several modes above 1300 cm \ By neglecting 


a number of terms in Eqn. (6.2), these modes can be removed from the spectral 


density with a negligible effect at smaller frequencies. Fig. 6.3 compares the dy¬ 


namics both with and without the high-frequency modes In both cases, the spectral 
densities have been scaled to give A = 35 cm -1 . 

On this short timescale, there are negligible differences between the predicted 


dynamics. Another interesting comparison is given by Fig. 6.4, which shows the 
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6.1. Structured Spectral Density 


dynamics using the Lorentz-Drude spectral density and the full structured spectral 


density, with a slightly greater discrepancy than Fig. 6.3[ though still a very small 


difference between the dynamics. 


(a) 




Figure 6.3: (a) Dynamics predicted using structured spectral density with (solid lines) and 
without (dashed lines) high-frequency modes, (b) Spectral densities used in the calculation 


of (a). 




















6.1. Structured Spectral Density 


(a) 




Figure 6.4: (a) Dynamics predicted using Lorentz-Drude (solid lines) and structured 
(dashed lines) spectral density, (b) Spectral densities used in the calculation of (a). 
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6.1. Structured Spectral Density 


It is possible that even a small difference in the short-time dynamics could lead 
to a larger difference in the longer-time dynamics, so it is important to follow the 
simulations to this limit. 

Fig. |6.5| repeats the two comparisons, but with long-time simulations. As before, 
there is excellent agreement for the structured spectral density and the same with 
high-frequency peaks omitted, and good (though less so than the former) agreement 
for the structured spectral density and the Lorentz-Drude. We now analyse these 
results. 


6.1.2 Discussion 


The effect of spectral structure on the dynamics of energy transfer in the FMO 
complex is very little at room temperature. 


We have seen in Fig. 6.3 that the high-frequency structure appears to be un¬ 


necessary in determining the efficiency of the energy transfer, and in Fig. 6.4 that 


the fully structured spectral density of Olbrich et al. 58 can be replaced with the 


Lorentz-Drude spectral density with only a small effect. 

This latter result is remarkable because it suggests that the actual character of 
the vibrational modes is unimportant, and that it is some other property of the bath 
that determines the energy transfer. It is interesting to look further into what this 
property is. 

Since we have shown that Forster theory provides a reasonable description of 
the dynamics, an investigation of absorption and fluorescence spectra, as in Section 


3.2.2, will be informative. 
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6.1. Structured Spectral Density 


(a) 



(b) 



Figure 6.5: Long-time dynamics for (a) fully structured (solid lines) vs. omitted high- 
frequency peaks (dashed lines) spectral densities, (b) Lorentz-Drude (solid lines) vs. fully 
structured (dashed lines) spectral densities. 


Fig. [6 )6] shows the fluorescence spectra of site 1 for both the Lorentz-Drude 


spectral density and the fully structured spectral density of Fig. 6.2 


The spectrum is narrowed to a small degree, but the effect of this on the spectral 
overlap is negligible compared to either of the following effects: the much greater 


narrowing that comes from decreasing the temperature (see Fig. 3.1) or the shift 


that comes about from increasing the reorganization energy. 
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6.1. Structured Spectral Density- 



Figure 6.6: Site 1 fluorescence spectra, F\[uj], for Lorentz-Drude (solid black line) and 
fully structured (dashed blue line) spectral densities. 


Fig. 6.7 shows the effect both on the dynamics and on F\ [u;] of increasing the 
reorganization energy from 35 cm -1 to 70 cm -1 , while using the same Lorentz-Drude 
spectral density: there is a more significant change in the dynamics, reflecting the 
fact that the change results in a larger effect on the position and width of the spectral 


line, and thus a larger effect on its overlap with other lines (and via Eqn. (3.25), on 
the energy transfer rate constants). 

These observations lead us to the conclusion that it is the temperature and the 
environmental reorganization energy that have the greatest effect on the dynamics 
of excitation energy transfer, rather than the precise shape of the spectral density 
JlV). 
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6.2. Static Disorder 


(a) 




Figure 6.7: (a) Energy transfer dynamics for the 7-site FMO complex using a Lorentz- 
Drude spectral density with A = 35 cm -1 (solid lines) and A = 70 cm -1 (dashed lines), 
(b) Site 1 fluorescence spectra, Fifo;], for the two reorganization energies. 


6.2 Static Disorder 

The phenomenon of static disorder causes inhomogeneous line-broadening in elec¬ 
tronic spectroscopy. That is, because each FMO complex in a sample is experiencing 
a different protein environment, giving different electronic transition energies, the 
spectral properties of each individual complex will differ and the superimposition of 
different spectral lines gives a broadened line. 

In modelling the effects of inhomogeneous broadening, we treat the site energies 
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6.2. Static Disorder 


as random variables with a Gaussian distribution 57 . 


Thus, from the original input data hies for our numerical simulations, a large 
number of other input hies are created, in which a random value is added to each 
site energy, drawn from a Gaussian distribution. The dynamics are found for each 
input hie, and the results are averaged. 

In order to carry out this procedure, we require a method of producing Gaussian 
random numbers. The standard Box-Muller algorithm was chosen for this purpose 
(Sec. 7.2 of |48j). 

The HEOM would be more difficult to use in finding the effect of this disorder, 
because a converged calculation for a single realization of Gaussian disorder would 
take a good fraction of an hour, whereas a Forster calculation is completed in a 
matter of seconds. We may require a large number of realizations to find the average 
effects, so it is preferable to use the faster technique, whose results have been shown 
to be accurate enough for this purpose at 300 K. 

The effect of inhomogeneous broadening on 2-dimensional electronic spectra of 
the FMO complex has been the subject of a number of studies (8,57, but in the 
following, we hope to find the effects, if any, on energy transfer dynamics in the 
hope once again of understanding the effect of the protein environment. 


6.2.1 Results 

An implementation of static disorder is characterized by the standard deviation, 
a, of the Gaussian distribution. Often, the full-width at half-maximum (FWHM), 
equal to 2\j2 In 2a, is specified instead. Several studies have used a FWHM on the 
order of 50 cm -1 (8,57,59 , and so this value will be used in this Section. 
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6.2. Static Disorder 


We show firstly the results for 3 realizations of static disorder in Fig. |6.8| along¬ 
side each set of results is a diagram showing the energy levels of the different sites 
for this realization, compared to the energy levels in the absence of disorder. 

Each individual set of site energies gives dynamics that differ a little from the 
case where disorder is neglected: this is a general observation that can be made 


for any one realization that is taken, although in some cases (as in Fig. 6.8 c), the 
differences are slightly more pronounced. 

In a real physical application, there will be a great number of FMO complexes, 
and the dynamics will be an average over many realizations. Fig. H compares 
the average dynamics for different numbers of realizations to the dynamics in the 
absence of disorder. A relatively small number of realizations is required for the 
average dynamics to converge. 


6.2.2 Discussion 

The results presented in this Section illustrate that inhomogeneous disordering ef¬ 
fects on the dynamics of the FMO complex are minimal. This can be explained by 


examining the energy level diagrams in the right-hand column of Fig. 6 


In the results we have observed, the populations of bacteriochlorophylls 5, 6 
and 7 remain close to zero, and only the other four have appreciable populations 
throughout the simulation. 

The electronic energy separations of the four sites involved in the transfer are 
larger than the differences induced by static disorder, so that the important transi¬ 
tion energies are not likely to be very different for two different realizations. 

Thus, while there will generally be small differences for the dynamics of these 
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6.2. Static Disorder 


Dynamics 
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Figure 6.8: Dynamics for 3 realizations of static disorder. The left-hand column gives the 
dynamics, and the right-hand column shows the site energies for these realizations. The 
solid line in each case is the situation in which static disorder is ignored, and the dashed 
line indicates its presence. 
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6.2. Static Disorder 


(a) 5 Realizations 



Figure 6.9: Average dynamics (dashed lines) for (a) 5 and (b) 50 realizations of static 
disorder. The solid lines in each case show the dynamics without disorder. 

two realizations, the most important features will be the same when comparing the 
two, and these small differences are not important in the FMO complex because of 
the effective averaging of the dynamics. 
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6.3. Conclusions 


6.3 Conclusions 

In the two subsections of this Chapter, conclusions have been drawn separately about 
the effects of a structured spectral density and of static disorder due to an inhomo¬ 
geneous environment. Here, we will make some comments about the implications 
for the FMO complex. 

Firstly, it was found that the structure of the spectral density is largely unimpor¬ 
tant at 300 K, and that the temperature and reorganization energy were far more 
important in determining the rate of energy transfer. 

Secondly, the effects of static disorder characterized by a physically reasonable 
FWHM were found to be very little, in addition to which the average behaviour 
converged with only a small number of realizations. 

This leads to a conclusion which may seem surprising at first: that is, the elec¬ 
tronic energy transfer dynamics in the FMO complex are largely insensitive to the 
fine-structure of the environment. Similarly, the disorder inherent in any experiment 
has little effect, giving a very robust transfer dynamics. 

The energy transfer dynamics of the FMO complex are very stable with respect 
to changes in its vibrational environment. Perhaps this is not a coincidence, but 
rather a result of the complex having evolved to perform its job so robustly. 
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Chapter 7 

Conclusions And Further Work 


In this project, we have investigated the physics of the Fenna-Matthews-Olson 
complex. The Hierarchical Equations of Motion were utilized to find numerically 
exact excitation energy transfer dynamics for the full 24-site trimer, with a fully 
converged bath correlation function. These are the first exact results that have been 
computed for this system. 

We have also applied the approximate Redfield and Fdrster methods (the latter 
of which neglects electronic coherence) to the problem, and found that at room 
temperature, Forster theory described the features of the dynamics well, leading to 
the conclusion that quantum-mechanical coherence is, for the most part, irrelevant 
for the efficiency of transfer. 

Having shown that Forster theory is applicable at this temperature, we then 
applied it to the problem of altering the effect of the environment on the system. 
It was found that the effects of adding structure to the spectral density of the 
bath phonons coupled to the electronic system or of including static disorder were 
minimal, and that the bath reorganization energies and the temperature were the 
most important parameters for determining the energy transfer dynamics. 

There has been much excitement recently about the importance of coherence in 
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7.1. Further Work 


photosynthetic systems (4, 36,60-62 , and the incoherent nature of Forster theory 
has led to its dismissal as an appropriate method for describing transfer (3 38 


However, the results of this work show that the search for an adequate description 
of electronic energy transfer need not necessarily require that the effects of coherence 
be included. 

Before moving on to discuss potential further work, it is quite satisfying to note 
that the energy transfer dynamics in the FMO complex, structurally characterized 
nearly 40 years ago |9], can be modelled in a very satisfactory way by a theory that 


is even more venerable 32 


7.1 Further Work 


The different facets of this work suggest a number of directions for further research: 

• The availability of numerically exact results for larger and more complex sys¬ 
tems than previously available allows effective benchmarking of approximate 
methods: the current glut of such methods, giving different results, means that 
it is a great advantage to be able to test them against a method (the HEOM) 
that can be used for many different systems and parameter regimes. 

• Having seen that we need not requite that coherence effects be included in 
an accurate prediction of energy transfer dynamics, methods can be sought 
that are more accurate than the Forster theory, and with greater ranges of 
applicability, while still being incoherent. 

• An approximate method that gives accurate results at lower temperatures 
might allow for a more comprehensive test of the effects of a structured spectral 
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7.1. Further Work 


density at these temperatures. 


With such investigations, perhaps we will be able to gain more insight into how 
important, if at all, quantum-mechanical interference effects are in biology. 
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Appendix A 


A.l The Bath Correlation Function 


The following correlation functions may be defined, with the tilde representing 
the interaction picture as in Chapter [2] and qj a being bath coordinates as in Section 

D 


(t) = (q ja (t)q ja (0))p , 

(A.la) 

CjJt) = ■ 

(A.lb) 

From this, using the fact that £j(f) = ^2 a gj a qja{t), and the independence of the 

bath degrees of freedom, it can be seen that: 


<xj(t) = (ut)u o)) = Y.^ c m 

a 

(A.2) 

The symmetrized and antisymmetrized correlation functions 

are given by (pp. 

122-124 of 28]): 


s Ja (t) = \ (c+(t) + Cjjt)), 

(A.3a) 

Ajjt) = 1 ( C%(t ) - crjt )). 

(A.3b) 
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A.l. The Bath Correlation Function 


From the definitions (A.la) and (A.lb): 


C ja {t) = tr B [q ja (0)q ja (t)e PHb ] 

Zb 

= -^tr B [q jQ {0)e-^ B e + ^ B q ja {t)e~^ B ] 
Zb 

Zb 

= -^-tr B [q ja (0)e~ mB q ja (t - i/3k)\ 

= (qja(t - iph)q ja (0))p 
= Cj~ a (t — i(3h). 


(A.4) 


Now, taking the Fourier transforms of the correlation functions, 


F[u] = / F(t)e lut dt. 


(A.5) 


and using Eqn. (A.4) gives: 


C7M = C- a {t)e^dt 


' —oo 

POO 


C+ a {t - if3h)e iujt dt 

C+ a {u)e iu]{ - u+iftn) du 

= I CUt)e^e-^dt 
J -oo 

= e-^C+M- 


' — OO 
POO 


' —oo 
POO 


(A.6) 


The Fourier transform of the antisymmetrized correlation function is then (using 
(A.3b) and (A.6)) given by Aj a [cu] = A (l — e - ^) C+ a [ui\, or: 


r + r _ ^Aj a [u] 
^ja Pi ^ _ e -pnu) ' 


(A.7) 
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A.l. The Bath Correlation Function 


Aj[u>\ can then be found using Eqn. (2.18) for qj a (t) and Eqn. (A.3b): 


A jaM = 2i I dt ' e Qja( 0)])^ 


' — OO 
roo 


— — dt ■ e 


ILUt 


2 i 


ik 




sin (ujjat) 


h 


AilJljaUJja ./ x 


dt ■ — e i ^ } ~ Uja ^ t ') 


irh 


‘2'lTTlrLCx^joi 


(5(u - LU ja ) - 5(u + Uja)) . 


Inserting this expression into (A.2) and using Eqns. (A.7) and (1.7) 


aj[uj\ = 


2vr 


_ £—/3Tiuj 

2nJj(uj) 

_ fD—fihu }' 


E 


h 


2 rrijaUja 


- LUja ) — 5(u + LUja)) 


(A.8) 


(A.9) 


Here, the inclusion of two delta functions means that the spectral density is 
defined for negative frequencies as well as positive, and that Jj(cu) = —Jj(—cu). By 
Fourier inversion: 

r °° Jj(Lu)e~ iu}t 


otj(t) — duu 


2 _ fD — pUu) 


(A.10) 


This expression can be rearranged to give Eqn. (2.27), making use of the fact 
that Jj(cu) is an odd function of c u. 
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A.2. The Pade Approximant 


A.2 The Pade Approximant 


In this appendix we find a rational function approximation to the Bose-Einstein 
function fsose^x) = (1 — e~ x )~ 1 . A similar procedure must be followed to find an 
approximation to the Fermi-Dirac function fFermiix) = (1 + e x ) _1 , and after the 
derivation is outlined, the results for both are given. Firstly, we take: 


1 1 e' 2 - e"*/ 2 1 e x / 2 + e~ x ' 2 1 1 , , . . 

I _ e -x - 2 e x / 2 - e~ x / 2 + 2e x / 2 - e ~ x / 2 _ 2 + 2 coth( A/ 2 )- 


(A.ll) 


Then, the continued-fraction representation of tanh(y) is given by the expression 


below 52 , where the notation is explained using the third expression: 


^ = A A ^ 


Inverting this expression gives coth(y): 


1 + JiA 
' 3 +.. 


.I / \ i . y y 2 y 2 

coth(y) = —|- 


y 3+5+7+...’ 


(A.12) 


(A.13) 


or, 


1 — e~ 


1 i V 4 x 7 4 x2 / a 

— _|_ l T . J— __'_ 

2 X 3+5+7+... 

1 i V 4 * 2 A ^A 

— ~ H-b x ■ -— -— --, 

2 x 6i+ b 2 + b 3 + ... 


(A.14) 


where b n — 2n + 1. The M th convergent of the continued fraction is 52 : 


C M (x ) = 


2n V 4 * 2 A * 2 A a m (x 2 ) 


&1+ &2 + 


■>M 


B m .(x 2 ) 


(A.15) 


A m(x 2 ) and Bm(x 2 ) are polynomials in x 2 . Each function, f M (x 2 ), is linked to 
/m- i(x 2 ) and to fM-2{x 2 ) by a recursion relation, and it is found [5l] that A 2 n{x 2 ) 
is of order N — 1, while B 2 n(x 2 ) is of order N. 
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A.2. The Pade Approximant 


Thus, C 2 n{x 2 ) is a rational function whose numerator has degree IV — 1 and 
whose denominator has degree N. It can be rewritten: 

nft 1 (* 2 + (j) N 


Con(x ) — 


n't. y 2 + sj) 


2r)jX 
x 2 + £ 2 ' 


(A.16) 


The roots of the numerator polynomial are then {—£ 2 } and those of the denom¬ 
inator are {—£ 2 }. The rjj coefficients (j = 1,2, ...,1V) are found by decomposing 


the first expression into partial fractions 52 


i i n Ar_1 (t 2 - b 2 i 

1 // 2 , t-2\r ( 2A 1 A7-J IIfc=l y^k S j 

Vj = o \( x +ij) C 2N{x ) x2= _ e = ^Nb N+ r 


(A.17) 


2 • — 2-— nhi(«-«l) 

Now, a straightforward method is required to find the roots of the two polyno¬ 


mials. It is shown in 52,63 that the convergent of a continued fraction can be 


expressed as, 


Cm{x 2 ) = [ £+ 2 ix = 


(A.18) 


11 


with the elements of D, D mn = b m 8 mn and of E, E mn = 8 m . n ±\, with each of m and 
n running from 1 to M. 

The roots of Bm(x 2 ) are the poles of Cm(x 2 ), and so are the roots of det (Z) + \ixE}j 

0. 

Thus, for B 2 n(x 2 ), the matrix A = D_~ l C_ D_ 1 is defined, with elements as in 


Eqn. (4.10), whose eigenvalues are {±2 /£l, ±2/£ 2 ,..., ±2/£jv} |52|. 


The roots of *4. 2 tv(^ 2 ) can be similarly found by using C 2 n(x 2 ) \ whose poles are 


these roots. A matrix A is defined, whose elements are also given in Eqn. (4.10), 
and whose eigenvalues are { 0 , ±2/( 1 , ±2/£ 2 , ..., ± 2 /£ jv - i }- 


Thus, the convergent C 2 w(^ 2 ) can be found using (A.16) and (A.17), so that the 
Pade approximant for the Bose-Einstein function is: 

fBose(x ) = \ + - + xC 2 n(x 2 ). (A. 19) 

2 x 














A.2. The Pade Approximant 


If a similar analysis is carried out for the Fermi-Dirac function, then the only differ¬ 
ence in defining the coefficient is that here b n — 2n — 1, and that: 


f Fermi (*^) ^ xC2n(% ) 


(A.20) 
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A.3. Supplementary Information 


A.3 Supplementary Information 


A.3.1 System Hamiltonians 


The system Hamiltonian for the 7-site FMO is that of Adolphs and Renger 13 
(energies in cm -1 ): 

( 


410 

-87.7 

5.5 

-5.9 

6.7 

-13.7 

-9.9 

-87.7 

530 

30.8 

8.2 

0.7 

11.8 

4.3 

5.5 

30.8 

210 

-53.5 

-2.2 

-9.6 

6.0 

-5.9 

8.2 

-53.5 

320 

-70.7 

-17.0 

-63.3 

6.7 

0.7 

-2.2 

-70.7 

480 

81.1 

-1.3 

-13.7 

11.8 

-9.6 

-17.0 

81.1 

630 

39.7 

-9.9 

4.3 

6.0 

-63.3 

-1.3 

39.7 

440 , 


(A.21) 


iltonian is given 

/ 

h A 

h 

h'\ 

=A 

=B 

=B 

h f 

h , 

h 

=B 

=A 

=B 

h 

h f 

V 

\=B 

=B 


(A.22) 


The matrix h A gives the intra-monomer couplings and site energies. The off- 









A.3. Supplementary Information 


diagonal elements (couplings) are (energies in cm x ): 


/ 

-80.3 

3.5 

-4.0 

4.5 

-10.2 

-4.9 

21.0 ^ 

-80.3 


23.5 

6.7 

0.5 

7.5 

1.5 

3.3 

3.5 

23.5 


-49.8 

-1.5 

-6.5 

1.2 

0.7 

-4.0 

6.7 

-49.8 


-63.4 

-13.3 

-42.2 

-1.2 

4.5 

0.5 

-1.5 

-63.4 


55.8 

4.7 

2.8 

-10.2 

7.5 

-6.5 

-13.3 

55.8 


33.0 

-7.3 

-4.9 

1.5 

1.2 

-42.2 

4.7 

33.0 


-8.7 

^ 21.0 

3.3 

0.7 

-1.2 

2.8 

-7.3 

-8.7 

/ 


The diagonal elements (site energies) found using the Olbrich et 
and Schmidt am Busch et al. (SaB, [l2|) studies are (energies in cm 


Site 

1 

2 

3 

4 

5 

6 

7 

8 

Olb 

186 

81 

0 

113 

65 

89 

492 

218 

SaB 

310 

230 

0 

180 

405 

320 

270 

505 


(A.23) 


l (Olb, (56)) 
): 

(A.24) 
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The matrix h gives the inter-monomer couplings (energies in cm 1 ): 


1.0 

0.3 

-0.6 

0.7 

2.3 

1.5 

0.9 

0.1 

1.5 

-0.4 

-2.5 

-1.5 

7.4 

5.2 

1.5 

0.7 

1.4 

0.1 

-2.7 

5.7 

4.6 

2.3 

4.0 

0.8 

0.3 

0.5 

0.7 

1.9 

-0.6 

-0.4 

1.9 

-0.8 

0.7 

0.9 

1.1 

-0.1 

1.8 

0.1 

-0.7 

1.3 

0.1 

0.7 

0.8 

1.4 

-1.4 

-1.5 

1.6 

-1.0 

0.3 

0.2 

-0.7 

4.8 

-1.6 

0.1 

5.7 

-2.3 

0.1 

0.6 

1.5 

-1.1 

4.0 

-3.1 

-5.2 

3.6 


(A.25) 
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A.3.2 Structured Spectral Density 


The structured spectral density in Fig. 6.1 is given by Eqn. (6.2), the parameters 


for which are given by the table below 58 


k 

Cjk/ cm -2 

ujjk/cmr 1 

7ifc/cm -1 

1 

7091.9081 

0.0000 

86.4772 

2 

2478.5629 

0.0000 

11.5643 

3 

617.8502 

1813.8341 

16.3525 

4 

2961.7657 

1752.8329 

22.4542 

5 

642.2013 

1692.3597 

15.1094 

6 

833.1906 

1600.5955 

21.9446 

7 

1108.6926 

1453.4383 

67.8879 

8 

4038.4676 

817.5591 

898.2804 

9 

957.3336 

1109.2920 

47.6728 

10 

683.7415 

721.8440 

51.0171 

11 

593.0216 

605.5993 

19.7612 

12 

1160.2597 

366.4331 

50.8997 

13 

97.4045 

972.4901 

17.8089 

14 

454.0769 

501.7510 

47.7413 

15 

725.2817 

243.7265 

59.1580 


(A.26) 


The spectral density of Fig. 6.3 (b) (dashed line) is given by dropping the last 


five entries in Eqn. 


(A.26) and scaling all rjjk values by a factor of 1.1194 (to give 


the correct reorganization energy). 
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